diff --git a/README.md b/README.md index f8f08d6..66aaa8e 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,6 @@

Logo

- # Porosity.jl [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ayushinav.github.io/Porosity.jl/stable/) diff --git a/docs/src/intro/getting_started.md b/docs/src/intro/getting_started.md index e18765a..298ddf8 100644 --- a/docs/src/intro/getting_started.md +++ b/docs/src/intro/getting_started.md @@ -93,8 +93,8 @@ resp = forward(m, []) println("size of T : ", size(T)) println("size of P : ", size(P)) -println("size of broadcasted array without f: ", size(T .+ P .+ dg .+ σ .+ ϕ .+ ρ .+ - T_solidus)) +println( + "size of broadcasted array without f: ", size(T .+ P .+ dg .+ σ .+ ϕ .+ ρ .+ T_solidus)) println("size of f : ", size(f)) println("size of Qinv (one of the outputs, averaged over frequency) : ", size(resp.Qinv)) ``` diff --git a/docs/src/intro/intro_figs.md b/docs/src/intro/intro_figs.md index 76fbc27..c71a2a0 100644 --- a/docs/src/intro/intro_figs.md +++ b/docs/src/intro/intro_figs.md @@ -77,7 +77,6 @@ See how the `plot_model` function returns a `figure` and an `axis`. If we want t arr3 = arr1 .^ 2 lines!(arr1, arr3) - ``` We only needed to pass `ax` to plot on the same `axis`. This might also seem a bit intuitive as well. At this point, we want to encourage creating empty `figure` beforehand and then always using mutating functions. This gives us more control on how the axes are arranged among other stuff. Here's an example to make things clear. diff --git a/docs/src/intro/intro_julia.md b/docs/src/intro/intro_julia.md index b603a35..5f17732 100644 --- a/docs/src/intro/intro_julia.md +++ b/docs/src/intro/intro_julia.md @@ -49,7 +49,7 @@ When functions and complex expressions are needed to be performed, `@.` can be u ```@repl x = [10, 20, 30, 40, 50] y = exp.(sin.(x ./ 10)) -y = @. exp(sin(x/10)) +y = @. exp(sin(x / 10)) ``` ## Functions @@ -70,8 +70,8 @@ Creating new functions is also similar to MATLAB, and is done via the following ```@repl function f(x, y) - a = x/(1 + x^2) - b = y/(1 + y^2) + a = x / (1 + x^2) + b = y / (1 + y^2) z = a + b return z end @@ -94,13 +94,13 @@ A lot of times, we want default arguments or named arguments. Consider the follo ```@example fn_demo function f(x, y, a=2; operation_type=2) if operation_type == 1 - return x+y/a + return x + y / a elseif operation_type == 2 - return x/(x^2 + a^2) + return x / (x^2 + a^2) elseif operation_type == 3 - return x/(x^2 + a^2) + return x / (x^2 + a^2) else - return y/a + return y / a end end ``` diff --git a/docs/src/solidus.md b/docs/src/solidus.md index f6c1203..2c45154 100644 --- a/docs/src/solidus.md +++ b/docs/src/solidus.md @@ -14,8 +14,8 @@ The dry solidus temperature can be obtained as a function of pressure `P` throug ```@example sol_plts f = Figure() -ax1 = Axis(f[1, 1]; xlabel="solidus temp (K)", ylabel="depth (km)", backgroundcolor=( - :magenta, 0.05)) +ax1 = Axis(f[1, 1]; xlabel="solidus temp (K)", + ylabel="depth (km)", backgroundcolor=(:magenta, 0.05)) ax1.yreversed = true xtfn(x) = @. string(round((x - 5) / 30.2; digits=2)) @@ -216,10 +216,14 @@ ps_nt = (; P, T, T_solidus, Ch2o, Cco2, D) size(ϕ) fig = Figure(; size=(800, 900)) -ax1 = Axis(fig[1, 1]; xlabel="bulk water conc. (ppm)", ylabel="bulk CO₂ conc. (ppm)", title="$(T[1]) K") -ax2 = Axis(fig[1, 2]; xlabel="bulk water conc. (ppm)", ylabel="bulk CO₂ conc. (ppm)", title="$(T[2]) K") -ax3 = Axis(fig[2, 1]; xlabel="bulk water conc. (ppm)", ylabel="bulk CO₂ conc. (ppm)", title="$(T[3]) K") -ax4 = Axis(fig[2, 2]; xlabel="bulk water conc. (ppm)", ylabel="bulk CO₂ conc. (ppm)", title="$(T[4]) K") +ax1 = Axis(fig[1, 1]; xlabel="bulk water conc. (ppm)", + ylabel="bulk CO₂ conc. (ppm)", title="$(T[1]) K") +ax2 = Axis(fig[1, 2]; xlabel="bulk water conc. (ppm)", + ylabel="bulk CO₂ conc. (ppm)", title="$(T[2]) K") +ax3 = Axis(fig[2, 1]; xlabel="bulk water conc. (ppm)", + ylabel="bulk CO₂ conc. (ppm)", title="$(T[3]) K") +ax4 = Axis(fig[2, 2]; xlabel="bulk water conc. (ppm)", + ylabel="bulk CO₂ conc. (ppm)", title="$(T[4]) K") crange = extrema(ϕ) cmap = :thermal @@ -232,7 +236,8 @@ contour!(ax3, vec(Ch2o), vec(Cco2), ϕ[3, :, :]; labels=true, color=:black, labe h = heatmap!(ax4, vec(Ch2o), vec(Cco2), ϕ[4, :, :]; colorrange=crange, colormap=cmap) contour!(ax4, vec(Ch2o), vec(Cco2), ϕ[4, :, :]; labels=true, color=:black, labelsize=14) Colorbar(fig[3, :], h; vertical=false, label="melt fraction") -Label(fig[0, :], "melt fraction at dry solidus temp : $(Int(round(T_solidus))) K"; fontsize=20) +Label(fig[0, :], "melt fraction at dry solidus temp : $(Int(round(T_solidus))) K"; + fontsize=20) nothing # hide ``` diff --git a/docs/src/tutorials/mixing_phases.md b/docs/src/tutorials/mixing_phases.md index b10dac4..5b382be 100644 --- a/docs/src/tutorials/mixing_phases.md +++ b/docs/src/tutorials/mixing_phases.md @@ -182,7 +182,6 @@ m = multi_phase_modelType(UHO2014, Dai_Karato2009, Sifre2014, HS_plus_multi_phas logsig_mat = zeros(length(T), length(ϕ1), length(ϕ2)) for i in eachindex(ϕ1), j in eachindex(ϕ2) - ϕ = [ϕ1[i], ϕ2[j]] ps_nt = (; ps_nt_..., ϕ) model = m(ps_nt) @@ -192,10 +191,14 @@ end fig = Figure(; size=(800, 900)) crange = extrema(logsig_mat) cmap = :thermal -ax1 = Axis(fig[1, 1]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[1]) K") -ax2 = Axis(fig[1, 2]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[2]) K") -ax3 = Axis(fig[2, 1]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[3]) K") -ax4 = Axis(fig[2, 2]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[4]) K") +ax1 = Axis(fig[1, 1]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[1]) K") +ax2 = Axis(fig[1, 2]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[2]) K") +ax3 = Axis(fig[2, 1]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[3]) K") +ax4 = Axis(fig[2, 2]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[4]) K") heatmap!(ax1, ϕ1, ϕ2, logsig_mat[1, :, :]; colorrange=crange, colormap=cmap) contour!(ax1, ϕ1, ϕ2, logsig_mat[1, :, :]; color=:black) @@ -229,7 +232,6 @@ m = multi_phase_modelType(UHO2014, Dai_Karato2009, Sifre2014, HS_minus_multi_pha logsig_mat = zeros(length(T), length(ϕ1), length(ϕ2)) for i in eachindex(ϕ1), j in eachindex(ϕ2) - ϕ = [ϕ1[i], ϕ2[j]] ps_nt = (; ps_nt_..., ϕ) model = m(ps_nt) @@ -239,10 +241,14 @@ end fig = Figure(; size=(800, 900)) crange = extrema(logsig_mat) cmap = :thermal -ax1 = Axis(fig[1, 1]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[1]) K") -ax2 = Axis(fig[1, 2]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[2]) K") -ax3 = Axis(fig[2, 1]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[3]) K") -ax4 = Axis(fig[2, 2]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[4]) K") +ax1 = Axis(fig[1, 1]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[1]) K") +ax2 = Axis(fig[1, 2]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[2]) K") +ax3 = Axis(fig[2, 1]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[3]) K") +ax4 = Axis(fig[2, 2]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[4]) K") heatmap!(ax1, ϕ1, ϕ2, logsig_mat[1, :, :]; colorrange=crange, colormap=cmap) contour!(ax1, ϕ1, ϕ2, logsig_mat[1, :, :]; color=:black) @@ -280,7 +286,6 @@ m = multi_phase_modelType(UHO2014, Dai_Karato2009, Sifre2014, GAL) logsig_mat = zeros(length(T), length(ϕ1), length(ϕ2)) for i in eachindex(ϕ1), j in eachindex(ϕ2) - ϕ = [ϕ1[i], ϕ2[j]] ps_nt = (; ps_nt_..., ϕ, m_GAL=[5.0, 4.0, 1.2]) model = m(ps_nt) @@ -290,10 +295,14 @@ end fig = Figure(; size=(800, 900)) crange = extrema(logsig_mat) cmap = :thermal -ax1 = Axis(fig[1, 1]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[1]) K") -ax2 = Axis(fig[1, 2]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[2]) K") -ax3 = Axis(fig[2, 1]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[3]) K") -ax4 = Axis(fig[2, 2]; xlabel="ol. vol. fraction", ylabel="opx. vol. fraction", title="Temp : $(T[4]) K") +ax1 = Axis(fig[1, 1]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[1]) K") +ax2 = Axis(fig[1, 2]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[2]) K") +ax3 = Axis(fig[2, 1]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[3]) K") +ax4 = Axis(fig[2, 2]; xlabel="ol. vol. fraction", + ylabel="opx. vol. fraction", title="Temp : $(T[4]) K") heatmap!(ax1, ϕ1, ϕ2, logsig_mat[1, :, :]; colorrange=crange, colormap=cmap) contour!(ax1, ϕ1, ϕ2, logsig_mat[1, :, :]; color=:black) diff --git a/docs/src/tutorials/stochastic_inverse.md b/docs/src/tutorials/stochastic_inverse.md index 2f6bf77..6d84745 100644 --- a/docs/src/tutorials/stochastic_inverse.md +++ b/docs/src/tutorials/stochastic_inverse.md @@ -78,7 +78,8 @@ f # hide Lets try to also infer the water content along with the temperature. Everything remains the same, except that now, in `ps_nt`, the vector corresponding to water content `Ch2o_ol` will be replaced by a corresponding distribution. ```@example rp_si -mdist = Poe2010Distribution(MvNormal([1200.0], [400.0;]), product_distribution([Uniform(50.0, 150.0)])) +mdist = Poe2010Distribution( + MvNormal([1200.0], [400.0;]), product_distribution([Uniform(50.0, 150.0)])) rdist = RockphyCondDistribution(Porosity.normal_dist) m_cache = mcmc_cache(mdist, rdist, 1000, Prior()); @@ -205,7 +206,8 @@ err_resp = multi_rp_response(RockphyCond(0.01 .* abs.(resp.cond.σ)), ps_nt_dist = (; T=product_distribution([Uniform(1200.0, 1400.0)]), P=[3.0], ρ=[3300.0], Ch2o_m=MvNormal([100.0], [20.0;]), ϕ=product_distribution([Uniform(0.01, 0.2)])) -m1 = multi_rp_modelDistributionType(SEO3Distribution, anharmonic_poroDistribution, Nothing, Nothing) +m1 = multi_rp_modelDistributionType( + SEO3Distribution, anharmonic_poroDistribution, Nothing, Nothing) mdist = m1(ps_nt_dist) # rdist = RockphyCondDistribution(Porosity.normal_dist) rdist = Porosity.multi_rp_responseDistribution(RockphyCondDistribution(normal_dist), @@ -283,8 +285,10 @@ ps_nt_dist = (; T=product_distribution([Uniform(1270.0, 1670.0)]), Ch2o_m=product_distribution([Uniform(50.0, 150.0)]), Cco2_m=[100.0], ϕ=product_distribution([Uniform(0.01, 0.2)]), P=[3.0], ρ=[3300.0]) -m_mixdist = two_phase_modelDistributionType(SEO3Distribution, Sifre2014Distribution, HS1962_plus) -m = multi_rp_modelDistributionType(typeof(m_mixdist), anharmonicDistribution, Nothing, Nothing) +m_mixdist = two_phase_modelDistributionType( + SEO3Distribution, Sifre2014Distribution, HS1962_plus) +m = multi_rp_modelDistributionType( + typeof(m_mixdist), anharmonicDistribution, Nothing, Nothing) mdist = m(ps_nt_dist); rdist = Porosity.multi_rp_responseDistribution(RockphyCondDistribution(normal_dist), diff --git a/src/models/anelastic/forward.jl b/src/models/anelastic/forward.jl index 7bdbf02..d4b9f5f 100644 --- a/src/models/anelastic/forward.jl +++ b/src/models/anelastic/forward.jl @@ -1,6 +1,5 @@ function SubsurfaceCore.forward(m::andrade_psp, p, params=default_params_andrade_psp) - @unpack n, β, τ_MR, E, G_UR, TR, PR, dR, Vstar, M, melt_alpha, ϕ_c, - elastic_type, params_elastic, melt_enhancement = params + @unpack n, β, τ_MR, E, G_UR, TR, PR, dR, Vstar, M, melt_alpha, ϕ_c, elastic_type, params_elastic, melt_enhancement = params resp_elastic = forward_for_anelastic(m, elastic_type, params_elastic) @@ -29,8 +28,7 @@ function SubsurfaceCore.forward(m::andrade_psp, p, params=default_params_andrade end function SubsurfaceCore.forward(m::eburgers_psp, p, params=default_params_eburgers_psp) - @unpack integration_params, elastic_type, params_elastic, params_btype, viscous_type, - params_viscous, JF10_visc, melt_enhancement = params + @unpack integration_params, elastic_type, params_elastic, params_btype, viscous_type, params_viscous, JF10_visc, melt_enhancement = params @unpack alf, DeltaB, DeltaP, sig = params_btype ω = 2.0f0π .* m.f @@ -40,18 +38,14 @@ function SubsurfaceCore.forward(m::eburgers_psp, p, params=default_params_eburge Ju = @. inv(G) - τ_maxwell, τ_L, - τ_H, - τ_P = calc_maxwell_times( + τ_maxwell, τ_L, τ_H, τ_P = calc_maxwell_times( G, m, params_btype, JF10_visc, params_viscous, viscous_type, melt_enhancement) J1_int_fn(x, ω) = x^(alf - 1) / (1 + (ω * x)^2) J2_int_fn(x, ω) = x^alf / (1 + (ω * x)^2) int1 = broadcast( - (l, - h, - omega) -> integrate_s(J1_int_fn, + (l, h, omega) -> integrate_s(J1_int_fn, omega, (l=l, h=h, N=integration_params.τ_integration_points, rule=integration_params.integration_method)), @@ -59,9 +53,7 @@ function SubsurfaceCore.forward(m::eburgers_psp, p, params=default_params_eburge τ_H, ω) int2 = broadcast( - (l, - h, - omega) -> integrate_s(J2_int_fn, + (l, h, omega) -> integrate_s(J2_int_fn, omega, (l=l, h=h, N=integration_params.τ_integration_points, rule=integration_params.integration_method)), @@ -79,8 +71,7 @@ function SubsurfaceCore.forward(m::eburgers_psp, p, params=default_params_eburge inv(x) * (exp(-0.5f0 * log(x / tau_p) * inv(sig))^2) * inv(1 + (ω * x)^2) end int11 = broadcast( - (omega, - tau_p) -> integrate_s((x, omega) -> J1_int_fn2(x, omega, tau_p), omega, + (omega, tau_p) -> integrate_s((x, omega) -> J1_int_fn2(x, omega, tau_p), omega, (l=eps(typeof(omega))^40, h=Inf, N=1, rule=Val(:quadgk)); rtol=1.0f16 * eps(eltype(ω))), ω, @@ -92,9 +83,8 @@ function SubsurfaceCore.forward(m::eburgers_psp, p, params=default_params_eburge (exp(-0.5f0 * log(x / tau_p) * inv(sig))^2) * inv(1 + (ω * x)^2) end int22 = broadcast( - (omega, - tau_p) -> integrate_s((x, omega) -> J2_int_fn2(x, omega, tau_p), omega, - (l=0.0f0, h=Inf, N=1, rule=Val(:quadgk))), + (omega, tau_p) -> integrate_s((x, omega) -> J2_int_fn2(x, omega, tau_p), + omega, (l=0.0f0, h=Inf, N=1, rule=Val(:quadgk))), ω, τ_P) # TODO : check this case @@ -114,7 +104,8 @@ function SubsurfaceCore.forward(m::eburgers_psp, p, params=default_params_eburge return RockphyAnelastic(J1, J2, Qinv, Ma, Va, Vave) end -function SubsurfaceCore.forward(m::premelt_anelastic, p, params=default_params_premelt_anelastic) +function SubsurfaceCore.forward( + m::premelt_anelastic, p, params=default_params_premelt_anelastic) @unpack params_xfit, elastic_type, params_elastic, viscous_params = params @unpack include_direct_melt_effect, β_B, poro_Λ, α_B, A_B, τ_pp = params_xfit @@ -163,8 +154,7 @@ function SubsurfaceCore.forward(m::premelt_anelastic, p, params=default_params_p end function SubsurfaceCore.forward(m::xfit_mxw, p, params=default_params_xfit_mxw) - @unpack α_a, α_b, α_c, α_τn, α2, β1, β2, τ_cutoff, melt_alpha, ϕ_c, - elastic_type, params_elastic, viscous_type, viscous_params = params + @unpack α_a, α_b, α_c, α_τn, α2, β1, β2, τ_cutoff, melt_alpha, ϕ_c, elastic_type, params_elastic, viscous_type, viscous_params = params resp_elastic = forward_for_anelastic(m, elastic_type, params_elastic) @unpack G, K, Vp, Vs = resp_elastic @@ -185,8 +175,8 @@ function SubsurfaceCore.forward(m::xfit_mxw, p, params=default_params_xfit_mxw) J_int_fn(x, _) = inv(x) * xfit_mxw_func(x, α_a, α_b, α_c, α2, β1, β2, α_τn, τ_cutoff) int1 = broadcast( - tau_norm_f -> integrate_s(J_int_fn, 0.0f0, ( - l=10.0f0^(-30.0f0), h=tau_norm_f, N=1, rule=Val(:quadgk))), + tau_norm_f -> integrate_s( + J_int_fn, 0.0f0, (l=10.0f0^(-30.0f0), h=tau_norm_f, N=1, rule=Val(:quadgk))), τ_norm_f) # TODO : check this case int2 = broadcast(J_int_fn, τ_norm_f, 0.0f0) @@ -204,9 +194,9 @@ function SubsurfaceCore.forward(m::xfit_mxw, p, params=default_params_xfit_mxw) return RockphyAnelastic(J1, J2, Qinv, Ma, Va, Vave) end -function SubsurfaceCore.forward(m::andrade_analytical, p, params=default_params_andrade_analytical) - @unpack α, β, η_ss, viscosity_method, viscosity_mech, elastic_type, - params_elastic, viscous_type, viscous_params = params +function SubsurfaceCore.forward( + m::andrade_analytical, p, params=default_params_andrade_analytical) + @unpack α, β, η_ss, viscosity_method, viscosity_mech, elastic_type, params_elastic, viscous_type, viscous_params = params resp_elastic = forward_for_anelastic(m, elastic_type, params_elastic) @unpack G, K, Vp, Vs = resp_elastic diff --git a/src/models/anelastic/utils.jl b/src/models/anelastic/utils.jl index 09ee9a9..2307388 100644 --- a/src/models/anelastic/utils.jl +++ b/src/models/anelastic/utils.jl @@ -1,8 +1,7 @@ # andrade_psp function calc_X̃(T, d, P, ϕ, params_anelastic) - @unpack n, β, τ_MR, E, G_UR, TR, PR, dR, Vstar, M, melt_alpha, ϕ_c, - elastic_type, params_elastic, melt_enhancement = params_anelastic + @unpack n, β, τ_MR, E, G_UR, TR, PR, dR, Vstar, M, melt_alpha, ϕ_c, elastic_type, params_elastic, melt_enhancement = params_anelastic X̃ = @. (d / dR)^(-M) * exp((-E / (gas_R * 1.0f3)) * (inv(T) - inv(TR)) - Vstar / (gas_R * 1.0f3) * (P / T - PR / TR) * 1.0f9) @@ -47,12 +46,7 @@ function get_η_diff(m, viscous_type::Val{HK2003}, params_viscous) x_ϕ_c_vec = get_melt_settings_for_x_ϕ_c(melt_enhancement) fH2O = @. calc_fH2O(m.Ch2o_ol, ch2o_o, P, m.T) ϵ_rate_diff = broadcast( - (T, - P, - σ, - d, - ϕ, - fH2O) -> sr_flow_law_calculation_HK2003( + (T, P, σ, d, ϕ, fH2O) -> sr_flow_law_calculation_HK2003( T, P * 1.0f9, σ, d, ϕ, fH2O, getfield(x_ϕ_c_vec, :diff), mechs, :diff), m.T, P, @@ -66,15 +60,15 @@ function get_η_diff(m, viscous_type::Val{HK2003}, params_viscous) end function get_η_diff(m, viscous_type::Val{xfit_premelt}, params_viscous) - resp_xfit_premelt = forward(xfit_premelt(m.T, m.P, m.dg, m.σ, m.ϕ, m.T_solidus), [], params_viscous) + resp_xfit_premelt = forward( + xfit_premelt(m.T, m.P, m.dg, m.σ, m.ϕ, m.T_solidus), [], params_viscous) return resp_xfit_premelt.η end function calc_maxwell_times(Gu, m::eburgers_psp, params_btype, JF10_visc, params_viscous, viscous_type, melt_enhancement) - @unpack TR, PR, dR, E, Vstar, Tau_LR, Tau_HR, Tau_MR, - Tau_PR, m_a, m_v, melt_alpha, ϕ_c = params_btype + @unpack TR, PR, dR, E, Vstar, Tau_LR, Tau_HR, Tau_MR, Tau_PR, m_a, m_v, melt_alpha, ϕ_c = params_btype x_ϕ_c = getfield(get_melt_settings_for_x_ϕ_c(melt_enhancement), :diff) if JF10_visc @@ -164,9 +158,7 @@ end # xfit_premelt aka premelt_anelastic function calc_Ap(Tn, ϕ, params) - @unpack α_B, - A_B, τ_pp, A_p_fac_1, A_p_fac_2, A_p_fac_3, σ_p_fac_1, σ_p_fac_2, σ_p_fac_3, - A_p_Tn_pts, σ_p_Tn_pts, include_direct_melt_effect, β, β_B, poro_Λ = params + @unpack α_B, A_B, τ_pp, A_p_fac_1, A_p_fac_2, A_p_fac_3, σ_p_fac_1, σ_p_fac_2, σ_p_fac_3, A_p_Tn_pts, σ_p_Tn_pts, include_direct_melt_effect, β, β_B, poro_Λ = params β_p = (include_direct_melt_effect) ? β : 0.0f0 @@ -189,9 +181,7 @@ function calc_Ap(Tn, ϕ, params) end function calc_σp(Tn, params) - @unpack α_B, - A_B, τ_pp, A_p_fac_1, A_p_fac_2, A_p_fac_3, σ_p_fac_1, σ_p_fac_2, σ_p_fac_3, - A_p_Tn_pts, σ_p_Tn_pts, include_direct_melt_effect, β, β_B, poro_Λ = params + @unpack α_B, A_B, τ_pp, A_p_fac_1, A_p_fac_2, A_p_fac_3, σ_p_fac_1, σ_p_fac_2, σ_p_fac_3, A_p_Tn_pts, σ_p_Tn_pts, include_direct_melt_effect, β, β_B, poro_Λ = params σ_p = 0.0f0 diff --git a/src/models/combine_models.jl b/src/models/combine_models.jl index 82654dc..5d1fff2 100644 --- a/src/models/combine_models.jl +++ b/src/models/combine_models.jl @@ -98,7 +98,8 @@ end function default_params(::Type{multi_rp_modelType{T1, T2, T3, T4}}) where {T1, T2, T3, T4} (; zip([:cond, :elastic, :visc, :anelastic], - [default_params(T1), default_params(T2), default_params(T3), default_params(T4)])...) + [default_params(T1), default_params(T2), + default_params(T3), default_params(T4)])...) end function SubsurfaceCore.from_nt(m::Type{T}, nt::NamedTuple) where {T <: multi_rp_modelType} diff --git a/src/models/conductivity/forward.jl b/src/models/conductivity/forward.jl index fc4a454..2845e5c 100644 --- a/src/models/conductivity/forward.jl +++ b/src/models/conductivity/forward.jl @@ -55,8 +55,7 @@ function SubsurfaceCore.forward(m::Jones2012, p, params=default_params_Jones2012 end function SubsurfaceCore.forward(m::Poe2010, p, params=default_params_Poe2010) - @unpack S_H100, H_H100, a_H100, r_H100, S_H010, H_H010, a_H010, r_H010, S_H001, H_H001, - a_H001, r_H001, S_A100, H_A100, S_A010, H_A010, S_A001, H_A001 = params + @unpack S_H100, H_H100, a_H100, r_H100, S_H010, H_H010, a_H010, r_H010, S_H001, H_H001, a_H001, r_H001, S_A100, H_A100, S_A010, H_A010, S_A001, H_A001 = params # Anhydrous σ_A100 = @. arrh_dry(S_A100, H_A100, boltz_k, m.T) diff --git a/src/models/elastic/forward.jl b/src/models/elastic/forward.jl index 967a7a4..c0c7d42 100644 --- a/src/models/elastic/forward.jl +++ b/src/models/elastic/forward.jl @@ -1,10 +1,9 @@ function SubsurfaceCore.forward(m::anharmonic, p, params=default_params_anharmonic) - @unpack T_K_ref, P_Pa_ref, Gu_0_ol, dG_dT, dG_dP, ν, Gu_0_crust, - dG_dT_crust, dG_dP_crust, Gu_tp_fn, Ku_tp_fn = params + @unpack T_K_ref, P_Pa_ref, Gu_0_ol, dG_dT, dG_dP, ν, Gu_0_crust, dG_dT_crust, dG_dP_crust, Gu_tp_fn, Ku_tp_fn = params - Gu₀, dG_dT₀, - dG_dP₀ = @. calc_Gu₀(Gu_0_ol, dG_dT, dG_dP, Gu_0_crust, dG_dT_crust, dG_dP_crust) #since χ is 1., we are always using ol + Gu₀, dG_dT₀, dG_dP₀ = @. calc_Gu₀( + Gu_0_ol, dG_dT, dG_dP, Gu_0_crust, dG_dT_crust, dG_dP_crust) #since χ is 1., we are always using ol ΔT = @. m.T - T_K_ref # K ΔP = @. m.P * 1.0f9 - P_Pa_ref # Pa @@ -19,7 +18,8 @@ function SubsurfaceCore.forward(m::anharmonic, p, params=default_params_anharmon return RockphyElastic(Gu_tp, Ku_tp, Vp, Vs) end -function SubsurfaceCore.forward(m::anharmonic_poro, p, params=default_params_anharmonic_poro) +function SubsurfaceCore.forward( + m::anharmonic_poro, p, params=default_params_anharmonic_poro) @unpack m_A, m_K, ν, p_anharmonic = params anh_p = forward(anharmonic(m.T, m.P, m.ρ), [], p_anharmonic) diff --git a/src/models/mixing_phases.jl b/src/models/mixing_phases.jl index 437a4e2..d002be2 100644 --- a/src/models/mixing_phases.jl +++ b/src/models/mixing_phases.jl @@ -99,7 +99,8 @@ function SubsurfaceCore.forward(model::two_phase_model{V, T1, T2, M}, @. σ1 = exp10(σ1) @. σ2 = exp10(σ2) - σ = broadcast((sig1, sig2, phi) -> mix_models([sig1, sig2], phi, model.mix), σ1, σ2, model.ϕ) + σ = broadcast( + (sig1, sig2, phi) -> mix_models([sig1, sig2], phi, model.mix), σ1, σ2, model.ϕ) return RockphyCond(log10.(σ)) end @@ -112,7 +113,8 @@ function SubsurfaceCore.forward(model::two_phase_model{V, T1, T2, M}, p, @. σ1 = exp10(σ1) @. σ2 = exp10(σ2) - σ = broadcast((sig1, sig2, phi) -> mix_models([sig1, sig2], phi, model.mix), σ1, σ2, model.ϕ) + σ = broadcast( + (sig1, sig2, phi) -> mix_models([sig1, sig2], phi, model.mix), σ1, σ2, model.ϕ) return RockphyCond(log10.(σ)) end @@ -359,7 +361,8 @@ function broadcast_helper_(m_tup, phi, mix, ::Val{2}) end function broadcast_helper_(m_tup, phi, mix, ::Val{3}) - broadcast((m1, m2, m3) -> mix_models((m1, m2, m3), phi, mix), m_tup[1], m_tup[2], m_tup[3]) + broadcast( + (m1, m2, m3) -> mix_models((m1, m2, m3), phi, mix), m_tup[1], m_tup[2], m_tup[3]) end function broadcast_helper_(m_tup, phi, mix, ::Val{4}) @@ -385,8 +388,8 @@ end function broadcast_helper_(m_tup, phi, mix, ::Val{8}) broadcast( - (m1, m2, m3, m4, m5, m6, m7, - m8) -> mix_models((m1, m2, m3, m4, m5, m6, m7, m8), phi, mix), + (m1, m2, m3, m4, m5, m6, m7, m8) -> mix_models( + (m1, m2, m3, m4, m5, m6, m7, m8), phi, mix), m_tup[1], m_tup[2], m_tup[3], @@ -401,13 +404,13 @@ function default_params(::Type{multi_phase_modelType{ T1, T2, T3, T4, T5, T6, T7, T8, M}}) where {T1, T2, T3, T4, T5, T6, T7, T8, M} (; zip([:m1, :m2, :m3, :m4, :m5, :m6, :m7, :m8], - [default_params(T1), default_params(T2), default_params(T3), - default_params(T4), default_params(T5), default_params(T6), - default_params(T7), default_params(T8)])...) + [default_params(T1), default_params(T2), default_params(T3), + default_params(T4), default_params(T5), default_params(T6), + default_params(T7), default_params(T8)])...) end -function SubsurfaceCore.from_nt(m::Type{T}, nt::NamedTuple) where {T <: - multi_phase_modelType} +function SubsurfaceCore.from_nt( + m::Type{T}, nt::NamedTuple) where {T <: multi_phase_modelType} ϕ = nt.ϕ m1 = m.types[1].parameters[1] m2 = m.types[2].parameters[1] @@ -429,8 +432,8 @@ function SubsurfaceCore.from_nt(m::Type{T}, nt::NamedTuple) where {T <: mix = from_nt(m.types[9].parameters[1], nt) - ϕ_vec = rearrange_ϕ(ϕ, multi_phase_modelType( - m1, m2, m3, m4, m5, m6, m7, m8, m.types[9].parameters[1])) + ϕ_vec = rearrange_ϕ( + ϕ, multi_phase_modelType(m1, m2, m3, m4, m5, m6, m7, m8, m.types[9].parameters[1])) return multi_phase_model( ϕ_vec, model1, model2, model3, model4, model5, model6, model7, model8, mix) diff --git a/src/models/pretty_printing.jl b/src/models/pretty_printing.jl index 90019cf..ee58e3a 100644 --- a/src/models/pretty_printing.jl +++ b/src/models/pretty_printing.jl @@ -8,8 +8,8 @@ const model_names_definition = ( Ch2o_cpx="Water concentration in clinopyroxene (ppm)", Cco2_m="CO₂ concentration in melt (ppm)", T_solidus="Solidus Temperature (K)") -function Base.show(io::IO, ::MIME"text/plain", m::model) where {model <: - AbstractRockphyModel} +function Base.show( + io::IO, ::MIME"text/plain", m::model) where {model <: AbstractRockphyModel} println(io, "Model : ", typeof(m).name.name) for k in propertynames(m) println(io, model_names_definition[k], " : ", getfield(m, k)) @@ -24,8 +24,8 @@ const resp_names_definition = ( J2="Imaginary part of dynamic compliance (Pa⁻¹)", Qinv="Attenuation", M="Modulus (Pa)", V="Anelastic S-wave velocity : (m/s)", Vave="Frequency averaged S-wave velocity (m/s)") -function Base.show(io::IO, ::MIME"text/plain", m::model) where {model <: - AbstractRockphyResponse} +function Base.show( + io::IO, ::MIME"text/plain", m::model) where {model <: AbstractRockphyResponse} println(io, "Rock physics Response : ", typeof(m).name.name) for k in propertynames(m) println(io, resp_names_definition[k], " : ", getfield(m, k)) diff --git a/src/models/thermals.jl b/src/models/thermals.jl index 425318a..35c0eba 100644 --- a/src/models/thermals.jl +++ b/src/models/thermals.jl @@ -112,7 +112,8 @@ function ΔT_h2o_Katz2003(ps_nt) etype_ = eltype(P .+ Ch2o_m) Ch2o_m_sat = @. etype_(12 * P^(0.6f0) + P) - Ch2o_m_ = @. ifelse(Ch2o_m * 1.0f-4 < Ch2o_m_sat, etype_(Ch2o_m * 1.0f-4), etype_(Ch2o_m_sat)) + Ch2o_m_ = @. ifelse( + Ch2o_m * 1.0f-4 < Ch2o_m_sat, etype_(Ch2o_m * 1.0f-4), etype_(Ch2o_m_sat)) dT = @. K * (Ch2o_m_)^γ T_solidus = @. T_solidus - dT return (; T_solidus) @@ -359,8 +360,10 @@ function f_melt(u, p, H2O_suppress_fn, CO2_suppress_fn) Ch2o_m = Ch2o * inv(D + u * (1 - D)) Cco2_m = get_Cco2_m((; ϕ=u, Cco2, Cco2_sat)).Cco2_m - T_new_H2O = H2O_suppress_fn((; T, T_solidus, Ch2o, Cco2, Cco2_sat, P, D, Ch2o_m)).T_solidus - T_new_CO2 = CO2_suppress_fn((; T, T_solidus, Ch2o, Cco2, Cco2_sat, P, D, Cco2_m)).T_solidus + T_new_H2O = H2O_suppress_fn((; + T, T_solidus, Ch2o, Cco2, Cco2_sat, P, D, Ch2o_m)).T_solidus + T_new_CO2 = CO2_suppress_fn((; + T, T_solidus, Ch2o, Cco2, Cco2_sat, P, D, Cco2_m)).T_solidus T_solidus_new = T_solidus - (2T_solidus - T_new_H2O - T_new_CO2) ΔT = max(zero(T - T_solidus_new), T - T_solidus_new) dTdF = -40 * P + 450 @@ -380,8 +383,8 @@ function get_melt_fraction_core( if f1 * f2 > 0 return zero(etype_) else - prob_init = IntervalNonlinearProblem{false}(f, ( - exp10(-15 * one(etype_)), one(etype_)), p) + prob_init = IntervalNonlinearProblem{false}( + f, (exp10(-15 * one(etype_)), one(etype_)), p) sol = solve(prob_init) return etype_(sol.u) end @@ -426,7 +429,8 @@ get_melt_fraction(ps_nt) (ϕ = [0.041345720244193085, 0.05481061353960702, 0.072300827427471, 0.09337629457529527, 0.1171677298066955],) ``` """ -function get_melt_fraction(ps_nt; H2O_suppress_fn=ΔT_h2o_Blatter2022, CO2_suppress_fn=ΔT_co2_Dasgupta2013) +function get_melt_fraction( + ps_nt; H2O_suppress_fn=ΔT_h2o_Blatter2022, CO2_suppress_fn=ΔT_co2_Dasgupta2013) ps_nt = (; Ch2o=0.0f0, Cco2=0.0f0, Cco2_sat=38.0f4, D=0.005f0, ps_nt...) @unpack Cco2, Cco2_sat, Ch2o, T, T_solidus, P, D = ps_nt diff --git a/src/models/viscous/forward.jl b/src/models/viscous/forward.jl index 73fe7b6..4b88245 100644 --- a/src/models/viscous/forward.jl +++ b/src/models/viscous/forward.jl @@ -29,12 +29,7 @@ function SubsurfaceCore.forward(m::HK2003, p, params=default_params_HK2003) for mech in keys(mechs) sr = broadcast( - (T, - P, - σ, - d, - ϕ, - fH2O) -> sr_flow_law_calculation_HK2003( + (T, P, σ, d, ϕ, fH2O) -> sr_flow_law_calculation_HK2003( T, P * 1.0f9, σ, d, ϕ, fH2O, getfield(x_ϕ_c_vec, mech), mechs, mech), m.T, P, diff --git a/src/probabilistic/models/combine_models.jl b/src/probabilistic/models/combine_models.jl index 65f5d58..dfeba10 100644 --- a/src/probabilistic/models/combine_models.jl +++ b/src/probabilistic/models/combine_models.jl @@ -24,7 +24,8 @@ Similar to `multi_rp_modelDistributionType` but instead accepts `Distribution` or `Nothing` ```julia -multi_rp_modelDistributionType()(SEO3Distribution, anharmonicDistribution, HK2003Distribution, Nothing) +multi_rp_modelDistributionType()( + SEO3Distribution, anharmonicDistribution, HK2003Distribution, Nothing) ``` Pass `Nothing` for the types you do not want responses of, eg. above @@ -60,7 +61,8 @@ Rock physics model to capture multiple rock physics model, susually constructed ## Usage ```julia -m = multi_rp_modelDistributionType()(SEO3Distribution, anharmonicDistribution, Nothing, Nothing) +m = multi_rp_modelDistributionType()( + SEO3Distribution, anharmonicDistribution, Nothing, Nothing) ps_nt_dist = (; T=product_distribution(Uniform(1200.0f0, 1400.0f0)), P=[3.0f0], ρ=product_distribution(Uniform(80.0f0, 120.0f0)), ϕ=[0.1f0]) model = m(ps_nt_dist) diff --git a/src/probabilistic/utils.jl b/src/probabilistic/utils.jl index 6e1cf0d..62bd0a6 100644 --- a/src/probabilistic/utils.jl +++ b/src/probabilistic/utils.jl @@ -138,7 +138,8 @@ function SubsurfaceCore.to_dist_nt(d::T) where {T <: multi_phase_modelDistributi end function SubsurfaceCore.to_dist_nt(d::T) where {T <: multi_rp_modelDistribution} - return merge(to_dist_nt(d.cond), to_dist_nt(d.elastic), to_dist_nt(d.visc), to_dist_nt(d.anelastic)) + return merge(to_dist_nt(d.cond), to_dist_nt(d.elastic), + to_dist_nt(d.visc), to_dist_nt(d.anelastic)) end function SubsurfaceCore.to_dist_nt(d::T) where {T <: multi_rp_responseDistribution} diff --git a/test/anelastic_test.jl b/test/anelastic_test.jl index 515af05..a15cd41 100644 --- a/test/anelastic_test.jl +++ b/test/anelastic_test.jl @@ -1,4 +1,4 @@ -@testitem "anelastic tests" tags = [:anelastic] begin +@testitem "anelastic tests" tags=[:anelastic] begin using JET T = collect(1073.0f0:30:1373.0f0) P = 2 .+ zero(T) @@ -97,7 +97,8 @@ @test_opt forward(model, []) @test_call forward(model, []) for k in fieldnames(RockphyAnelastic) - @test all(isapprox.(log10.(getfield(out_, k)), log10.(getfield(outs[m], k)), rtol=1.5f-2),) + @test all(isapprox.( + log10.(getfield(out_, k)), log10.(getfield(outs[m], k)), rtol=1.5f-2),) end modelD = methodsD_list[i](inps[m]...) @test sample_type(modelD) <: methods_list[i] diff --git a/test/combine_models_test.jl b/test/combine_models_test.jl index cfabe79..03c130c 100644 --- a/test/combine_models_test.jl +++ b/test/combine_models_test.jl @@ -20,7 +20,8 @@ @test resp.elastic.Vp == resp_anharmonic.Vp @test resp.elastic.Vs == resp_anharmonic.Vs - m_dist = multi_rp_modelDistributionType(SEO3Distribution, anharmonicDistribution, Nothing, Nothing) + m_dist = multi_rp_modelDistributionType( + SEO3Distribution, anharmonicDistribution, Nothing, Nothing) modelD = m_dist(ps_nt) sample_type(modelD) <: multi_rp_modelType diff --git a/test/cond_test.jl b/test/cond_test.jl index e39d578..cf9f060 100644 --- a/test/cond_test.jl +++ b/test/cond_test.jl @@ -1,4 +1,4 @@ -@testitem "conductivity tests" tags = [:cond] begin +@testitem "conductivity tests" tags=[:cond] begin using JET T = collect(1273.0f0:30:1573.0f0) diff --git a/test/elastic_test.jl b/test/elastic_test.jl index 17f3048..4b63659 100644 --- a/test/elastic_test.jl +++ b/test/elastic_test.jl @@ -1,4 +1,4 @@ -@testitem "elastic tests" tags = [:elastic] begin +@testitem "elastic tests" tags=[:elastic] begin using JET T = collect(1273.0f0:30:1573.0f0) ρ = collect(3300.0f0:100.0f0:4300.0f0) diff --git a/test/phase_mixing_test.jl b/test/phase_mixing_test.jl index f288b1f..f39b4eb 100644 --- a/test/phase_mixing_test.jl +++ b/test/phase_mixing_test.jl @@ -1,4 +1,4 @@ -@testitem "two_phase" tags = [:phase_mixing] begin +@testitem "two_phase" tags=[:phase_mixing] begin using JET m1 = two_phase_modelType(SEO3, Gaillard2008, HS1962_plus) ps_nt = (; T=[1200.0f0, 1400.0f0] .+ 273, P=3.0f0, ρ=3300.0f0, Ch2o_m=100.0f0, ϕ=0.1f0) @@ -13,7 +13,8 @@ resp = forward(model, []) @test :σ ∈ propertynames(resp) - m1dist = two_phase_modelDistributionType(SEO3Distribution, Gaillard2008Distribution, HS1962_plus) + m1dist = two_phase_modelDistributionType( + SEO3Distribution, Gaillard2008Distribution, HS1962_plus) modelD = m1dist(ps_nt) @test sample_type(modelD) <: two_phase_modelType @@ -24,7 +25,7 @@ @test_call forward_helper(sample_type(modelD), ps_nt, [], (; σ=no_tf), params_) end -@testitem "multi_phase" tags = [:phase_mixing] begin +@testitem "multi_phase" tags=[:phase_mixing] begin using JET m1 = multi_phase_modelType(SEO3, Sifre2014, Zhang2012, HS_minus_multi_phase) ps_nt = (; T=[1200.0f0] .+ 273, Ch2o_ol=[1.0f0], Ch2o_m=[1000.0f0], diff --git a/test/tune_rp_test.jl b/test/tune_rp_test.jl index 89f2519..766f922 100644 --- a/test/tune_rp_test.jl +++ b/test/tune_rp_test.jl @@ -12,7 +12,8 @@ @report_opt m(ps_nt) @report_call m(ps_nt) - m_dist_ = two_phase_modelDistributionType(SEO3Distribution, Gaillard2008Distribution, HS1962_plus) + m_dist_ = two_phase_modelDistributionType( + SEO3Distribution, Gaillard2008Distribution, HS1962_plus) m_dist = tune_rp_modelDistributionType(fn_list, m_dist_) modelD = m_dist(ps_nt) diff --git a/test/visc_test.jl b/test/visc_test.jl index 7a0db01..bc3fa31 100644 --- a/test/visc_test.jl +++ b/test/visc_test.jl @@ -1,4 +1,4 @@ -@testitem "viscosity tests" tags = [:visc] begin +@testitem "viscosity tests" tags=[:visc] begin using JET T = collect(1073.0f0:30:1373.0f0) P = 2 .+ zero(T) @@ -15,11 +15,14 @@ outs = ( HZK2011=RockphyViscous( - Float32[1.0f-12, 2.0f-12, 5.0f-12, 1.2f-11, 2.8f-11, 6.2f-11, 1.36f-10, 2.86f-10, 5.86f-10, 1.1710001f-9, 2.283f-9], + Float32[1.0f-12, 2.0f-12, 5.0f-12, 1.2f-11, 2.8f-11, 6.2f-11, + 1.36f-10, 2.86f-10, 5.86f-10, 1.1710001f-9, 2.283f-9], [8.9619f18, 3.8154f18, 1.6603f18, 7.410f17, 3.397f17, 1.601f17, 7.75f16, 3.85f16, 1.96f16, 1.03f16, 5.5f15]), HK2003=RockphyViscous( - Float32[3.0f-11, 7.999999f-11, 1.8999999f-10, 4.5999998f-10, 1.05f-9, 2.3499998f-9, 5.09f-9, 1.0709999f-8, 2.191f-8, 4.3609997f-8, 8.4579995f-8], + Float32[ + 3.0f-11, 7.999999f-11, 1.8999999f-10, 4.5999998f-10, 1.05f-9, 2.3499998f-9, + 5.09f-9, 1.0709999f-8, 2.191f-8, 4.3609997f-8, 8.4579995f-8], [2.3787f17, 1.0128f17, 4.408f16, 1.968f16, 9.03f15, 4.26f15, 2.06f15, 1.03f15, 5.2f14, 2.8f14, 1.5f14]), xfit_premelt=RockphyViscous(zero(T), @@ -34,7 +37,8 @@ @test_opt forward(model, []) @test_call forward(model, []) for k in fieldnames(RockphyViscous) - @test all(isapprox.(log10.(getfield(out_, k)), log10.(getfield(outs[m], k)), rtol=1.0f-2),) + @test all(isapprox.( + log10.(getfield(out_, k)), log10.(getfield(outs[m], k)), rtol=1.0f-2),) end modelD = methodsD_list[i](inps[m]...) @test sample_type(modelD) <: methods_list[i]