Merge pull request #156 from nichollsh/hn/morels

Hn/morels

Author: nichollsh

.zenodo.json

@@ -22,7 +22,7 @@
         }
     ],
     "title": "AGNI: An open-source model for extreme atmospheres on rocky exoplanets",
-    "version": "1.8.0",
+    "version": "1.8.1",
     "access_right": "open",
     "related_identifiers": [
         {

CITATION.cff

@@ -8,9 +8,9 @@ authors:
 
 title: "AGNI: An open-source model for extreme atmospheres on rocky exoplanets"
 doi: 10.21105/joss.07726
-date-released: 2025-11-03
+date-released: 2025-11-20
 
-version: 1.8.0
+version: 1.8.1
 url: "https://github.com/nichollsh/AGNI"
 type: software
 license: "GPL-3.0"

Project.toml

@@ -1,7 +1,7 @@
 name = "AGNI"
 uuid = "ede838c1-9ec3-4ebe-8ae8-da4091b3f21c"
 authors = ["Harrison Nicholls <harrison.nicholls@physics.ox.ac.uk>"]
-version = "1.8.0"
+version = "1.8.1"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"

codemeta.json

@@ -13,11 +13,11 @@
   "identifier": "",
   "codeRepository": "https://github.com/nichollsh/AGNI",
   "datePublished": "2024-12-29",
-  "dateModified": "2025-10-21",
+  "dateModified": "2025-11-20",
   "dateCreated": "2024-09-13",
   "description": "An open-source model for extreme atmospheres on rocky exoplanets.",
   "keywords": "exoplanets, radiation, astronomy, atmosphere, optimization, julia, fortran, convection, chemistry, magma, planetary, physics",
   "license": "GPL-3.0",
   "title": "AGNI",
-  "version": "1.8.0"
+  "version": "1.8.1"
 }

misc/utilities/grid_agni.jl

@@ -34,18 +34,17 @@ mass_arr::Array{Float64, 1} = 10.0 .^ vcat( range(start=log10(0.5),  stop=log10(
 #    enter the least-important parameters first
 grid::OrderedDict = OrderedDict{String,Array{Float64,1}}((
 
+    "frac_core"     =>       range(start=0.2,   stop=0.7,   step=0.1),
+    "frac_atm"      =>       range(start=0.00,  stop=0.15,  step=0.03),
+    "mass_tot"      =>       mass_arr,  # M_earth
+
     # metallicities here are by MASS fraction relative to hydrogen (converted to mole below)
-    "metal_C"       => 10 .^ range(start=-1.0,  stop=3.0,   step=2.0),
     # "metal_S"       => 10 .^ range(start=-1.0,  stop=3.0,     step=2.0),
     # "metal_O"       => 10 .^ range(start=-1.0,  stop=3.0,     step=2.0),
-
-    "frac_core"     =>       range(start=0.2,   stop=0.7,   step=0.1),
-    "frac_atm"      =>       range(start=0.00,  stop=0.15,  step=0.03),
+    "metal_C"       => 10 .^ range(start=-1.0,  stop=3.0,   step=2.0),
 
     # "Teff"          =>       range(start=2500,  stop=6000,  step=700.0),
     "instellation"  => 10 .^ range(start=log10(1.0),  stop=log10(2500.0),  length=5), # S_earth
-
-    "mass_tot"      =>       mass_arr,  # M_earth
 ))
 
 # Variables to record
@@ -56,17 +55,16 @@ output_keys =  ["succ", "flux_loss",
 
 # Grid management options
 save_netcdfs = false        # NetCDF file for each case
-save_plots   = false        # plots for each case
-save_ncdf_tp = true         # a single NetCDF containing all T(p) solutions
+save_plots   = true         # plots for each case
 
 # Runtime options
-AGNI.solver.ls_increase= 1.02
-modwrite::Int          = 1           # frequency to write CSV file
-modplot::Int           = 2            # Plot during runtime (debug)
+AGNI.solver.ls_increase= 1.1
+modwrite::Int          = 2            # frequency to write CSV file
+modplot::Int           = 0            # Plot during runtime (debug)
 frac_min::Float64      = 0.001        # 0.001 -> 1170 bar for Earth
 frac_max::Float64      = 1.0
 transspec_p::Float64   = 2e3    # Pa
-fc_floor::Float64      = 500.0   # K
+fc_floor::Float64      = 300.0   # K
 
 # =============================================================================
 # Parse keys and flatten grid
@@ -84,6 +82,18 @@ mkdir(output_dir)
 # Backup config to output dir
 cp(joinpath(ROOT_DIR,cfg_base), joinpath(output_dir,"base_config.toml"))
 
+# IO directory
+# IO_DIR::String = "/tmp/"
+# if haskey(ENV,"TMPDIR")
+#     IO_DIR = abspath(ENV["TMPDIR"])
+# end
+# IO_DIR = joinpath(IO_DIR,"agni_grid") * "/"
+IO_DIR = output_dir
+@info "IO folder: $IO_DIR"
+
+mkdir(joinpath(output_dir,"nc"))
+mkdir(joinpath(output_dir,"pl"))
+
 # Results path
 result_table_path::String = joinpath(output_dir,"result_table.csv")
 result_emits_path::String = joinpath(output_dir,"result_emits.csv")
@@ -340,6 +350,7 @@ atmosphere.setup!(atmos, ROOT_DIR, output_dir,
                                 p_top,
                                 mf_dict, "";
 
+                                IO_DIR=IO_DIR,
                                 condensates=condensates,
                                 κ_grey_lw=grey_lw,
                                 κ_grey_sw=grey_sw,
@@ -418,7 +429,8 @@ function update_structure!(atmos, mass_tot, frac_atm, frac_core)
     atmos.grav_surf = phys.grav_accel(atmos.interior_mass, atmos.rp)
 
     # surface pressure
-    atmos.p_boa = mass_atm * atmos.grav_surf / (4 * pi * atmos.rp^2)
+    atmos.p_oboa = mass_atm * atmos.grav_surf / (4 * pi * atmos.rp^2)
+    atmos.p_boa = atmos.p_oboa
     atmosphere.generate_pgrid!(atmos)
 end
 
@@ -523,16 +535,14 @@ for (i,p) in enumerate(grid_flat)
 
     # Set temperature array based on interpolation from last solution
     max_steps = Int(cfg["execution"]["max_steps"])
-    easy_start = Bool(cfg["execution"]["easy_start"])
-    if i > 1
-        if succ_last
-            # last iter was successful
-            setpt.fromarrays!(atmos, result_profs[i-1]["p"], result_profs[i-1]["t"])
-        else
-            # last iter failed -> restore initial guess for T(p)
-            setpt.request!(atmos, cfg["execution"]["initial_state"])
-            easy_start = true
-        end
+    if succ_last && (i > 1)
+        # last iter was successful
+        setpt.fromarrays!(atmos, result_profs[i-1]["p"], result_profs[i-1]["t"])
+        easy_start = false
+    else
+        # last iter failed -> restore initial guess for T(p)
+        setpt.request!(atmos, cfg["execution"]["initial_state"])
+        easy_start = Bool(cfg["execution"]["easy_start"])
     end
 
     # Solve for RCE
@@ -554,20 +564,22 @@ for (i,p) in enumerate(grid_flat)
                                             save_frames=false,
                                             radiative_Kzz=false,
                                             perturb_all=perturb_all,
-
                                             )
 
     # Report radius
     @info @sprintf("    found r_phot = %.3f R_earth",atmos.transspec_r/R_earth)
 
     # Write NetCDF file for this case
     if save_netcdfs
-        save.write_ncdf(atmos, joinpath(atmos.OUT_DIR,@sprintf("%08d.nc",i)))
+        save.write_ncdf(atmos, joinpath(atmos.OUT_DIR,"nc",@sprintf("%08d.nc",i)))
     end
 
     # Make plot for this case
     if save_plots
-        plotting.plot_pt(atmos, joinpath(atmos.OUT_DIR,@sprintf("%08d_pt.png",i)))
+        plotting.combined(plotting.plot_pt(atmos,""), plotting.plot_fluxes(atmos, ""),
+                            plotting.plot_vmr(atmos,""), plotting.plot_radius(atmos, ""),
+                            "Index = $i     Success = $succ",
+                            joinpath(atmos.OUT_DIR,"pl",@sprintf("%08d_pl.png",i)))
     end
 
     # Record keys (all in SI)

res/config/structure_grid.toml

@@ -24,21 +24,21 @@ title = "Base config for structure calcs"
 [composition]
     p_surf          = 100.0
     p_top           = 1e-5
-    metallicities   = {"C"=0.01, "S"=0.1, "O"=0.1}  # mass fraction
+    metallicities   = {"C"=0.1, "O"=10.0}  # mass fraction
     condensates     = []
 
 [execution]
     clean_output    = false
     verbosity       = 1
     max_steps       = 50
     max_runtime     = 200
-    num_levels      = 35
+    num_levels      = 30
     solution_type   = 3
     solver          = "newton"
-    dx_max          = 100.0
-    initial_state   = ["loglin","800"]
+    dx_max          = 300.0
+    initial_state   = ["iso","2000"]
     linesearch      = 1
-    easy_start      = true
+    easy_start      = false
     grey_start      = false
     converge_atol   = 1.0
     converge_rtol   = 1.0e-2

src/atmosphere.jl

@@ -33,7 +33,7 @@ module atmosphere
     import ..spectrum
 
     # Code versions
-    const AGNI_VERSION::String     = "1.8.0"  # current agni version
+    const AGNI_VERSION::String     = "1.8.1"  # current agni version
     const SOCVER_minimum::Float64  = 2407.2    # minimum required socrates version
 
     # Hydrostatic and gravity calc (constants and limits)

src/solver.jl

@@ -117,7 +117,7 @@ module solver
     plateau_r::Float64 =    0.98    # Cost ratio for determining whether to increment plateau_i
     #    linesearch
     ls_tau::Float64    =    0.7     # backtracking downscale size
-    ls_increase::Float64 =  1.08    # factor by which the cost can increase from last step before triggering linesearch
+    ls_increase::Float64 =  0.9    # factor by which the cost can increase from last step before triggering linesearch
     ls_max_steps::Int    =  10      # maximum steps undertaken by linesearch routine
     ls_min_scale::Float64 = 1.0e-5  # minimum step scale allowed by linesearch
     ls_max_scale::Float64 = 0.99    # maximum step scale allowed by linesearch
@@ -323,14 +323,15 @@ module solver
             # Set new temperatures
             _set_tmps!(x)
 
-            # Do saturation aloft here, only
-            # chemistry.calc_composition!(atmos, oceans, chem, rainout)
+            # Do saturation aloft here, only. Keep chemistry fixed.
             if rainout
                 # reset back to post-chemistry mixing ratios
                 for g in atmos.gas_names
                     @. atmos.gas_vmr[g] = atmos.gas_cvmr[g]
                 end
                 chemistry._sat_aloft!(atmos)
+            else
+                atmosphere.calc_layer_props!(atmos)
             end
 
             # Calculate fluxes
@@ -383,7 +384,6 @@ module solver
             ok::Bool = true
 
             # Evalulate residuals at x
-            # chemistry.calc_composition!(atmos, oceans, chem, rainout)
             ok = ok && _fev!(x, resid)
 
             # For each level...
@@ -779,7 +779,7 @@ module solver
                             end
 
                             ls_cost = _ls_func(ls_alpha)
-                            if ls_cost <= c_cur*ls_increase
+                            if ls_cost <= c_cur#*ls_increase
                                 # this scale is good enough
                                 break
                             end

test/runtests.jl

@@ -740,7 +740,7 @@ if suite > 15
 
     # temperature profile
     @info "Temperature..."
-    arr_e[:] .= [160.66339503503497, 172.09710160055477, 178.03603947440752, 179.20888073828587, 179.91013709468695, 179.28194461691936, 179.48335170955022, 181.60786079397423, 183.11232837720024, 182.74336048754452, 183.44034094824806, 186.23392659907887, 190.33872200527122, 195.3550015009265, 201.1036011339729, 207.48982773345082, 214.28148826113636, 221.70556437388842, 231.68725754907447, 246.4815220524279, 263.35634452191607, 287.5627049103804, 320.30148032756995, 354.8531587447293, 392.11917744416724, 429.26323831006886, 460.58872296988125, 487.0625306433956, 514.4554623098902, 543.1684695093304, 570.3876923217258, 593.8467013136508, 610.9314215719224, 620.7713434612795, 625.5983452566984, 628.122351599638, 629.6697885689832, 630.5708014622505, 630.9679678191192, 631.0915922829595, 631.121625484581, 631.1282688098997, 631.1296751975232, 631.129812085828, 631.1298132423226]
+    arr_e[:] .= [160.663417882603, 172.097128812409, 178.036068253899, 179.208908569656, 179.910164720311, 179.281972493295, 179.483379941445, 181.607889256078, 183.112356682626, 182.74338943682, 183.440372062087, 186.233959168495, 190.338756315127, 195.355039534191, 201.103646189501, 207.489881120782, 214.2815536312, 221.705657433829, 231.687457669314, 246.481687372921, 263.356595002906, 287.5630763156, 320.301515497962, 354.852774535303, 392.118464544988, 429.26259724383, 460.588392363893, 487.062314373208, 514.455301328919, 543.168345970076, 570.387595711342, 593.846621505299, 610.931353333751, 620.771283792757, 625.59828953617, 628.122297051458, 629.66973541692, 630.570749507789, 630.967916271267, 631.091540930136, 631.121574633071, 631.128221749457, 631.12965996539, 631.129771581873, 631.129771953692]
     if all(abs.(arr_e[:] .- arr_o_tmp[:]) .< abs.(arr_e[:].*rtol) .+ 0.5)
         @info "Pass"
     else
@@ -773,7 +773,7 @@ if suite > 15
 
     # eddy diffusion coefficients
     @info "Kzz..."
-    arr_e[:] .= [297631.1773891829, 245779.83944484807, 202961.69913189596, 167603.05242102614, 138404.3556049965, 114292.46289807544, 94381.18488545873, 77938.71821912714, 64360.74949697297, 53148.24480646298, 43889.108627309775, 36243.03799898127, 29929.015295112327, 24714.979924151103, 20409.299358103035, 16853.726022315273, 13917.581189404838, 11492.952116772693, 9490.725907098406, 7837.314323464027, 6471.949185555334, 5344.448944073614, 13566.607073176909, 12453.628747107967, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415, 8814.614880555415]
+    arr_e[:] .= [293515.330261667, 242381.028037358, 200155.006214063, 165285.323017843, 136490.40572033, 112711.948729345, 93076.0174630662, 76860.9284503441, 63470.7251477871, 52413.2746248124, 43282.1800995569, 35741.8445533194, 29515.1364634408, 24373.2043307382, 20127.0656526946, 16620.6612101877, 13725.1193905086, 11334.0197421416, 9359.48168174701, 7728.93460077965, 6382.45066279981, 5270.54226322015, 13528.8039526634, 12294.5144847605, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378, 8732.90905977378]
     if all(abs.(arr_e[:] .- arr_o_Kzz[:]) .< abs.(arr_e[:].*rtol) .+ 1e4)
         @info "Pass"
     else