Merge pull request #243 from brown-ccv/241-accept-vector-for-chamber-…

…arguments

241 accept vector for chamber arguments

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.8.5"
 manifest_format = "2.0"
-project_hash = "5417bfc0838607d8e63d898ada8e96159e5679e5"
+project_hash = "307cac97db0b7035e7716e25208fc63e328a2438"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
@@ -808,6 +808,12 @@ git-tree-sha1 = "c13304c81eec1ed3af7fc20e75fb6b26092a1102"
 uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 version = "0.3.2"
 
+[[deps.Memoize]]
+deps = ["MacroTools"]
+git-tree-sha1 = "2b1dfcba103de714d31c033b5dacc2e4a12c7caa"
+uuid = "c03570c3-d221-55d1-a50c-7939bbd78826"
+version = "0.4.4"
+
 [[deps.Missings]]
 deps = ["DataAPI"]
 git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272"

Project.toml

@@ -12,6 +12,7 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 IterableTables = "1c8ee90f-4401-5389-894e-7a04a3dc0f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"

scripts/runcode.jl

@@ -65,7 +65,7 @@ julia> output_dirname = "MyDirname"
 julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth, output_dirname)
 ```
 """
-function chamber(
+@memoize function chamber(
     composition::Union{Silicic,Mafic},
     end_time::Float64,
     log_volume_km3::Float64,
@@ -100,7 +100,7 @@ function chamber(
     )
 
     to = TimerOutput()
-    @timeit to "chamber" begin
+    @timeit to "vol$(log_volume_km3)_h2o$(InitialConc_H2O)_gas$(InitialConc_CO2)_vfr$(log_vfr)_dep$(depth)" begin
         rc = rheol_composition_dict[composition_str]
         param = Param{Float64}(;
             composition=composition,
@@ -248,10 +248,126 @@ function chamber(
             dtmax=odesetting.max_step,
         )
     end
-    @info(to)
+    print_timer_log(io, to)
     close(io)
     df = DataFrame(sol)
     write_csv(df, erupt_saved, path)
     plot_figs(df, path)
     return df
 end
+
+"""
+    chamber(composition::Union{Silicic,Mafic}, end_time::Float64, log_volume_km3_vector::Union{Float64,Vector{Float64}}, InitialConc_H2O_vector::Union{Float64,Vector{Float64}}, InitialConc_CO2_vector::Union{Float64,Vector{Float64}}, log_vfr_vector::Union{Float64,Vector{Float64}}, depth_vector::Union{Float64,Vector{Float64}}, output_dirname::String)
+
+Simulate the eruption of a volcano using a model for the frequency of eruptions of upper crustal magma chambers based on Degruyter and Huber (2014).
+
+## Arguments
+- `composition`: The magma composition. Use `Silicic()` for rhyolite composition (arc magma) or `Mafic()` for basalt composition (rift).
+- `end_time`: Maximum magma chamber evolution duration in seconds.
+- `log_volume_km3`: The initial volume of the chamber in logarithmic scale. The actual initial chamber volume is calculated as 10^(`log_volume_km3`) in km³.
+- `InitialConc_H2O`: The initial weight fraction of water in the magma (exsolved + dissolved).
+- `InitialConc_CO2`: The initial weight fraction of CO₂ in the magma (exsolved + dissolved).
+- `log_vfr`: Magma recharge rate in km³/yr calculated as 10^(`log_vfr`).
+- `depth`: Depth of the magma chamber in meters.
+- `output_dirname`(optional): Name of the output directory. Defaults to current timestamp.
+
+## Returns
+A `DataFrame` containing the solution with columns:
+- `time`: Simulation timestamps in seconds.
+- `P+dP`: Pressure in Pa.
+- `T`: Temperature in K.
+- `eps_g`: Gas volume fraction.
+- `V`: Volume of the magma chamber in m³.
+- `rho_m`: Density of the melt in kg/m³.
+- `rho_x`: Density of magma crystal in kg/m³.
+- `X_CO2`: Mole fraction of CO2 in the gas.
+- `total_mass`: Total mass of magma chamber in kg.
+- `total_mass_H2O`: Total mass of water in the magma in kg.
+- `total_mass_CO2`: Total mass of CO₂ in the magma in kg.
+
+## Outputs
+A directory named after `output_dirname` or the default value, containing the following files:
+- `out.csv`: a CSV file containing the solution columns listed above.
+- `eruptions.csv`, A CSV file containing the datas of eruptions with the following columns: time_of_eruption (sec), duration_of_eruption (sec), mass_erupted (kg) and volume_erupted (km³).
+- Figures for P+dP(t), T(t), eps_g(t), V(t), X_CO2(t), total_mass(t).
+
+## References
+- W. Degruyter and C. Huber (2014). A model for eruption frequency of upper crustal silicic magma chambers. Earth Planet. Sci. Lett. https://doi.org/10.1016/j.epsl.2014.06.047
+
+## Examples
+```
+# Run a simulation with silicic magma chamber
+julia> composition = Silicic()
+julia> end_time = 3e9
+julia> log_volume_km3 = 0.2
+julia> InitialConc_H2O = [0.04, 0.045]
+julia> InitialConc_CO2 = [0.001, 0.0012]
+julia> log_vfr = -3.3
+julia> depth = 8e3
+
+julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth)
+
+# Run a simulation with mafic magma chamber, with custom directory name "MyDirname"
+julia> composition = Mafic()
+julia> end_time = 3e9
+julia> log_volume_km3 = 0.2
+julia> InitialConc_H2O = [0.01, 0.012]
+julia> InitialConc_CO2 = 0.001
+julia> log_vfr = -3.3
+julia> depth = [7e3, 8e3]
+julia> output_dirname = "MyDirname"
+
+julia> chamber(composition, end_time, log_volume_km3, InitialConc_H2O, InitialConc_CO2, log_vfr, depth, output_dirname)
+```
+"""
+function chamber(
+    composition::Union{Silicic,Mafic},
+    end_time::Float64,
+    log_volume_km3_vector::Union{Float64,Vector{Float64}},
+    InitialConc_H2O_vector::Union{Float64,Vector{Float64}},
+    InitialConc_CO2_vector::Union{Float64,Vector{Float64}},
+    log_vfr_vector::Union{Float64,Vector{Float64}},
+    depth_vector::Union{Float64,Vector{Float64}},
+    output_dirname::String=get_timestamp(),
+)::String
+    path0 = joinpath(pwd(), output_dirname)
+    mkdir(path0)
+    to = TimerOutput()
+    @timeit to "chamber" begin
+        for log_volume_km3 in log_volume_km3_vector,
+            InitialConc_H2O in InitialConc_H2O_vector,
+            InitialConc_CO2 in InitialConc_CO2_vector,
+            log_vfr in log_vfr_vector,
+            depth in depth_vector
+
+            dataset = "vol$(log_volume_km3)_h2o$(InitialConc_H2O)_gas$(InitialConc_CO2)_vfr$(log_vfr)_dep$(depth)"
+            chamber(
+                composition,
+                end_time,
+                log_volume_km3,
+                InitialConc_H2O,
+                InitialConc_CO2,
+                log_vfr,
+                depth,
+                joinpath(output_dirname, dataset),
+            )
+        end
+    end
+    io0 = open("$path0/$output_dirname.log", "w+")
+    with_logger(SimpleLogger(io0)) do
+        @info(
+            "Arguments:",
+            composition,
+            end_time,
+            log_volume_km3_vector,
+            InitialConc_H2O_vector,
+            InitialConc_CO2_vector,
+            log_vfr_vector,
+            depth_vector,
+            path0
+        )
+    end
+    print_timer_log(io0, to)
+    close(io0)
+    return output_dirname
+end

src/Chamber.jl

@@ -4,6 +4,7 @@ using DataFrames
 using Dates
 using DifferentialEquations
 using LinearAlgebra
+using Memoize
 using Parameters
 using Plots
 using Roots
@@ -24,7 +25,8 @@ include("./ic_finder-utils.jl")
 include("./ic_finder.jl")
 include("./utils-matrix.jl")
 include("../scripts/runcode.jl")
-export Silicic,
+export @memoize,
+    Silicic,
     Mafic,
     QNDF,
     FBDF,
@@ -36,8 +38,11 @@ export Silicic,
     TimerOutput,
     @timeit,
     reset_timer!,
+    print_timer,
+    print_timer_log,
     SimpleLogger,
     global_logger,
+    with_logger,
     get_timer,
     SciMLBase,
     VectorContinuousCallback,

src/utils.jl

@@ -4,7 +4,7 @@ struct Silicic end
 struct Mafic end
 
 function get_timestamp()::String
-    return Dates.format(now(), "YYmmddHHMM")
+    return Dates.format(now(), "YYmmddHHMMSSsss")
 end
 
 """
@@ -400,6 +400,11 @@ function record_erupt_end(
     return nothing
 end
 
+function print_timer_log(io::IO, to::TimerOutput)::Nothing
+    print_timer(io, to; compact=true, allocations=false, sortby=:firstexec)
+    return nothing
+end
+
 """
     ic_phase_conversion(phase_here::T, composition::Union{Silicic,Mafic}, M_h2o::T, M_co2::T, M_tot::T, P::T, Temp::T, V::T, rho_m::T, param_IC::ParamICFinder{T})::NamedTuple{(:eps_g_temp, :X_co2_temp, :C_co2_temp, :phase),NTuple{4,T}} where {T<:Float64}