Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
<p align="center">
<img src="docs/src/assets/logo.png" alt="Logo" width="20%">
</p>

# Porosity.jl

[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ayushinav.github.io/Porosity.jl/stable/)
Expand Down
4 changes: 2 additions & 2 deletions docs/src/intro/getting_started.md
Original file line number Diff line number Diff line change
Expand Up @@ -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))
```
Expand Down
1 change: 0 additions & 1 deletion docs/src/intro/intro_figs.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
14 changes: 7 additions & 7 deletions docs/src/intro/intro_julia.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
```
Expand Down
19 changes: 12 additions & 7 deletions docs/src/solidus.md
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand All @@ -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
```

Expand Down
39 changes: 24 additions & 15 deletions docs/src/tutorials/mixing_phases.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down
12 changes: 8 additions & 4 deletions docs/src/tutorials/stochastic_inverse.md
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down Expand Up @@ -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),
Expand Down Expand Up @@ -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),
Expand Down
42 changes: 16 additions & 26 deletions src/models/anelastic/forward.jl
Original file line number Diff line number Diff line change
@@ -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)

Expand Down Expand Up @@ -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
Expand All @@ -40,28 +38,22 @@ 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)),
τ_L,
τ_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)),
Expand All @@ -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(ω))),
ω,
Expand All @@ -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

Expand All @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down
24 changes: 7 additions & 17 deletions src/models/anelastic/utils.jl
Original file line number Diff line number Diff line change
@@ -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)
Expand Down Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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

Expand Down
Loading