Skip to content

Commit

Permalink
Merge pull request #68 from nichollsh/eos
Browse files Browse the repository at this point in the history
Groundwork for reading AQUA equation of state
  • Loading branch information
nichollsh authored Sep 19, 2024
2 parents fc7a084 + 43a50fe commit bcd6343
Show file tree
Hide file tree
Showing 8 changed files with 882 additions and 573 deletions.
1 change: 1 addition & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ PCHIPInterpolation = "afe20452-48d1-4729-9a8b-50fb251f06cd"
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
ScatteredInterpolation = "3f865c0f-6dca-5f4d-999b-29fe1e7e3c92"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

Expand Down
4 changes: 3 additions & 1 deletion docs/paper/paper.md
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,9 @@ bibliography: paper.bib

It is important that we are able to accurately model the atmospheres of (exo)planets. This is because their atmospheres play a crucial role in setting the environment and conditions on the planet, including how the planet evolves over astronomical timescales. Additionally, it is primarily by observation of their atmospheres that we are able to characterise exoplanets. There is demand for accurate atmosphere models in the context of lava worlds: planets with permanent or fleeting magma oceans.

AGNI[^1] is a Julia program designed to solve for the temperature and radiation environment within the atmospheres of rocky (exo)planets. It leverages a well established FORTRAN code to calculate radiative fluxes from a given atmospheric temperature structure and composition, which -- alongside representations of convection and other processes -- enables an energy-conserving numerical solution for the atmospheric conditions. In contrast to most other numerical atmosphere codes, AGNI uses an Newton-Raphson optimisation method to obtain this solution which enables improved performance and scalability. AGNI was specifically developed for use alongside planetary interior models within a coupled framework, although it can also be easily applied to scientific problems as a standalone code. It can be interacted with as a library, or as an executable which reads configuration files. There are notebook tutorials in the repository, and the documentation can be read online [here](https://nichollsh.github.io/AGNI/).
AGNI[^1] is a Julia program designed to solve for the temperature and radiation environment within the atmospheres of rocky (exo)planets. It leverages a well established FORTRAN code to calculate radiative fluxes from a given atmospheric temperature structure and composition, which -- alongside representations of convection and other processes -- enables an energy-conserving numerical solution for the atmospheric conditions. In contrast to most other numerical atmosphere codes, AGNI uses an Newton-Raphson optimisation method to obtain this solution which enables improved performance and scalability. AGNI was specifically developed for use alongside planetary interior models within a coupled framework, although it can also be easily applied to scientific problems as a standalone code.

AGNI can be interacted with as a library; it is used in this way within the [Jupyter notebook tutorials](https://github.com/nichollsh/AGNI/tree/main/tutorials) in the repository. It can also be used as an executable program, where it reads TOML configuration files from the disk and outputs figures and NetCDF data to a specified directory. The documentation can be read online [here](https://nichollsh.github.io/AGNI/).

[^1]: AGNI can be found on GitHub [here](https://github.com/nichollsh/AGNI).

Expand Down
1,250 changes: 704 additions & 546 deletions misc/blame_emission.ipynb

Large diffs are not rendered by default.

42 changes: 21 additions & 21 deletions misc/rce-runaway.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@
"nlev_centre = 35\n",
"p_surf = 270.0 # bar\n",
"t_surf = 3000.0\n",
"p_top = 1e-5 # bar \n",
"p_top = 1e-5 # bar\n",
"mole_fractions = Dict([(\"H2O\", 1.0), (\"N2\",1.0e-9)]) # needs a small amount of background gas with zero opacity\n",
"condensates = [\"H2O\"]\n",
"\n",
Expand Down Expand Up @@ -133,10 +133,10 @@
"source": [
"# Setup atmosphere\n",
"atmos = atmosphere.Atmos_t()\n",
"atmosphere.setup!(atmos, ROOT_DIR, output_dir, \n",
"atmosphere.setup!(atmos, ROOT_DIR, output_dir,\n",
" spectral_file,\n",
" instellation, 0.375, 0.0, 48.19,\n",
" t_surf, \n",
" t_surf,\n",
" gravity, radius,\n",
" nlev_centre, p_surf, p_top,\n",
" mole_fractions, \"\",\n",
Expand Down Expand Up @@ -243,14 +243,14 @@
}
],
"source": [
"solver_success = nl.solve_energy!(atmos, \n",
"solver_success = nl.solve_energy!(atmos,\n",
" sol_type=1, # Conserve energy, but with fixed surface temperature\n",
" sens_heat=true, # Do not include sensible heat transport\n",
" latent=true, # Include condensation\n",
" method=1, # Use the Newton-Raphson method\n",
" dx_max=400.0, # Allow large step sizes because of the poor initial guess\n",
" ls_method=1 , # Enable Linesearch\n",
" save_frames=false, modplot=1, # disable plotting \n",
" save_frames=false, modplot=1, # disable plotting\n",
" conv_atol=0.1\n",
" )\n",
"println(\"Solver success? $solver_success\")\n",
Expand Down Expand Up @@ -489,32 +489,32 @@
"atmos = deepcopy(base_atm)\n",
"\n",
"for (i,tmp) in enumerate(tmp_arr)\n",
" # Set new surface temperature \n",
" # Set new surface temperature\n",
" @printf(\"Running model for tmp_surf = %.1f K \\n\",tmp)\n",
" atmos.tmp_surf = tmp \n",
" atmos.tmp_surf = tmp\n",
"\n",
" # Help solver by changing initial guess\n",
" clamp!(atmos.tmpl, 0.0, tmp)\n",
" clamp!(atmos.tmp, 0.0, tmp)\n",
"\n",
" # Run model\n",
" solver_success = nl.solve_energy!(atmos, \n",
" sol_type=1, \n",
" sens_heat=true, \n",
" latent=true, \n",
" method=1, \n",
" solver_success = nl.solve_energy!(atmos,\n",
" sol_type=1,\n",
" sens_heat=true,\n",
" latent=true,\n",
" method=1,\n",
" dx_max=100.0, # Smaller steps\n",
" ls_method=1 , \n",
" ls_method=1 ,\n",
" save_frames=false, modplot=1,\n",
" modprint=0,\n",
" conv_atol=0.1 \n",
" conv_atol=0.1\n",
" )\n",
"\n",
" # Store result\n",
" push!(atm_arr, deepcopy(atmos))\n",
"\n",
" @printf(\"--------------------------------- \\n\")\n",
"end \n",
"end\n",
"println(\"Done!\")"
]
},
Expand Down Expand Up @@ -545,15 +545,15 @@
"ylims = (arr_P[1], arr_P[end])\n",
"yticks = 10.0 .^ round.(Int,range( log10(ylims[1]), stop=log10(ylims[2]), step=1))\n",
"\n",
"plt = plot(framestyle=:box, size=(700,500), dpi=300, \n",
"plt = plot(framestyle=:box, size=(700,500), dpi=300,\n",
" leg=:topright, legcolumn=-1,\n",
" tickfontsize=fontsize, guidefontsize=fontsize, \n",
" tickfontsize=fontsize, guidefontsize=fontsize,\n",
" legendfontsize=fontsize, leg_title=\"Ts [K]\")\n",
"\n",
"p = Plots.palette(:viridis, num_loops)\n",
"for (i,this_atm) in enumerate(atm_arr)\n",
" plot!(plt, this_atm.tmpl, this_atm.pl* 1.0e-5, lc=p[i], linewidth=linewidth, label=@sprintf(\"%d\",this_atm.tmp_surf))\n",
"end \n",
"end\n",
"\n",
"xlabel!(plt, \"Temperature [K]\")\n",
"xaxis!(plt, minorgrid=true)\n",
Expand Down Expand Up @@ -583,7 +583,7 @@
"for this_atm in atm_arr\n",
" push!(olr_arr, this_atm.flux_u_lw[1])\n",
" push!(tot_arr, this_atm.flux_tot[1])\n",
"end \n"
"end\n"
]
},
{
Expand Down Expand Up @@ -617,15 +617,15 @@
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.10.4",
"display_name": "Julia 1.10.5",
"language": "julia",
"name": "julia-1.10"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.4"
"version": "1.10.5"
}
},
"nbformat": 4,
Expand Down
63 changes: 63 additions & 0 deletions res/config/L-98-59d.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
# AGNI configuration file
title = "L 98-59 d"

[planet]
tmp_surf = 2700.0
instellation = 81690.0
albedo_b = 0.0
s0_fact = 0.25
zenith_angle = 54.74
surface_material= "res/surface_albedos/lunar_marebasalt.dat"
radius = 6.59051e6
gravity = 16.73
tmp_int = 0.0
turb_coeff = 0.001
wind_speed = 2.0

[files]
input_sf = "res/spectral_files/Honeyside/256/Honeyside.sf"
input_star = "res/stellar_spectra/l-98-59.txt"
output_dir = "out/"

[composition]
p_top = 1e-6
p_surf = 215.82518354535077e3
vmr_dict = {H2S= 0.23988329190194904, SO2= 0.0044668359215096305, N2= 5.623413251903491e-6, H2= 0.7093907341517219}
include_all = true
chemistry = 1
condensates = []

[execution]
clean_output = true
verbosity = 1
max_steps = 20000
max_runtime = 400
num_levels = 60
continua = true
rayleigh = true
cloud = false
aerosol = false
overlap_method = 4
thermo_funct = true
sensible_heat = false
latent_heat = true
convection = true
solution_type = 3
solvers = ["newton"]
dx_max = 400.0
initial_state = ["loglin", "700"]
linesearch = 0
easy_start = false
converge_atol = 1.0e-2
converge_rtol = 1.0e-3

[plots]
at_runtime = true
temperature = true
fluxes = true
contribution = true
emission = true
albedo = true
mixing_ratios = true
animate = true

2 changes: 2 additions & 0 deletions src/AGNI.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ module AGNI
include("plotting.jl")
include("energy.jl")
include("solver.jl")
include("realgas.jl")

# Import src files
import .phys
Expand All @@ -31,6 +32,7 @@ module AGNI
import .plotting
import .energy
import .solver
import .realgas

# Export
# export atmosphere
Expand Down
10 changes: 5 additions & 5 deletions src/atmosphere.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1326,12 +1326,12 @@ module atmosphere
# evaluate gamma at band centre, converting from m to nm
ga = _gamma(1.0e9 * atmos.bands_cen[i])

# calculate dh reflectance (eq 3 from Hammond24)
atmos.surf_r_arr[i] = (1-ga) / (1+2*ga*mu)

# calculate emissivity (eq 4 from Hammond24 and eq 15.29 from Hapke 2012)
atmos.surf_e_arr[i] = 1.0 - ((1-ga)/(1+ga)) * (1- ga/(3*(1+ga)))
# calculate spherical reflectance (eq 4 from Hammond24)
atmos.surf_r_arr[i] = (1-ga)/(1+ga) * (1- ga/(3*(1+ga)))
end

# calculate emissivity (eq 4 from Hammond24 and eq 15.29 from Hapke 2012)
@turbo @. atmos.surf_e_arr = 1.0 - atmos.surf_r_arr
end

#######################################
Expand Down
83 changes: 83 additions & 0 deletions src/realgas.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
# Contains routines for loading and working with real-gas equations of state

# Not for direct execution
if (abspath(PROGRAM_FILE) == @__FILE__)
thisfile = @__FILE__
error("The file '$thisfile' is not for direct execution")
end

module realgas

import ..phys

using ScatteredInterpolation # 2D interpolation
using LoggingExtras
using DelimitedFiles

mutable struct Aqua_t

# Axis
npts_t::Int
npts_p::Int

# Array limits
lim_prs::Array{Float64}
lim_tmp::Array{Float64}

# Interpolators
itp_rho # Interpolator for density [kg/m^3]
itp_gad # Interpolator for adiabatic gradient (dlog(T)/dlog(P))_S
itp_mmw # Interpolator for mmw [kg/mol]

# Struct
Aqua_t() = new()
end

function read_aqua(fpath::String)::Aqua_t

# Read the file
@debug "Reading AQUA table at '$fpath'"
data::Array{Float64,2} = readdlm(fpath, Float64; header=false, skipstart=19)

# Parse the columns
# input
arr_prs::Array{Float64} = data[:,1]
arr_tmp::Array{Float64} = data[:,2]
# output
arr_rho::Array{Float64} = data[:,3]
arr_gad::Array{Float64} = data[:,4]
arr_mmw::Array{Float64} = data[:,8]

# Create struct
aqua = Aqua_t()

# Dimensions are hardcoded inside the file
aqua.npts_p = 1093
aqua.npts_t = 301

# Check dimensions vs data
if length(arr_prs) != aqua.npts_p*aqua.npts_t
@error "AQUA lookup table has invalid dimensions"
end

# Limits
aqua.lim_prs = [minimum(arr_prs), maximum(arr_prs)]
aqua.lim_tmp = [minimum(arr_tmp), maximum(arr_tmp)]

# Create interpolators
eval_pts = vcat(arr_prs', arr_tmp')
aqua.itp_rho = interpolate(NearestNeighbor(), eval_pts, arr_rho)
aqua.itp_gad = interpolate(NearestNeighbor(), eval_pts, arr_gad)
aqua.itp_mmw = interpolate(NearestNeighbor(), eval_pts, arr_mmw)

return aqua
end

# Evaluate density from temperature and pressure
function eval_aqua_rho(aqua::Aqua_t, prs::Float64, tmp::Float64)::Float64
prs = max(min(prs, aqua.lim_prs[2]), aqua.lim_prs[1])
tmp = max(min(tmp, aqua.lim_tmp[2]), aqua.lim_tmp[1])
return aqua.itp_rho(tmp,prs)
end

end

0 comments on commit bcd6343

Please sign in to comment.