Skip to content

Commit 15d471e

Browse files
author
Andrew Cheng-An Hsieh
committed
Implement deep atmospheric heating module
1 parent f560f1f commit 15d471e

9 files changed

Lines changed: 1295 additions & 12 deletions

File tree

‎misc/test_deep_heating.ipynb‎

Lines changed: 343 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,343 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"id": "4f297f36",
7+
"metadata": {},
8+
"outputs": [],
9+
"source": [
10+
"# Setup environment\n",
11+
"using Pkg\n",
12+
"Pkg.activate(\"..\")\n",
13+
"\n",
14+
"ENV[\"RAD_DIR\"] = abspath(\"..\", \"socrates\")\n",
15+
"\n",
16+
"using Revise\n",
17+
"include(\"../src/AGNI.jl\")\n",
18+
"\n",
19+
"using Printf\n",
20+
"using Plots"
21+
]
22+
},
23+
{
24+
"cell_type": "code",
25+
"execution_count": null,
26+
"id": "0b394f17",
27+
"metadata": {},
28+
"outputs": [],
29+
"source": [
30+
"# Common parameters\n",
31+
"ROOT_DIR = abspath(\"..\")\n",
32+
"OUT_DIR_BASE = abspath(\"..\", \"out\", \"test_deep_heating\")\n",
33+
"spectral_file = \"res/spectral_files/Dayspring/48/Dayspring.sf\"\n",
34+
"star_file = \"res/stellar_spectra/sun.txt\"\n",
35+
"\n",
36+
"# Planet parameters (hot Jupiter-like)\n",
37+
"instellation = 5000.0 # W/m² (moderate stellar flux)\n",
38+
"tmp_surf = 1500.0 # K\n",
39+
"gravity = 10.0 # m/s²\n",
40+
"radius = 7.0e7 # m (roughly Jupiter radius)\n",
41+
"p_surf = 100.0 # bar\n",
42+
"p_top = 1e-5 # bar\n",
43+
"nlev = 100 # number of levels\n",
44+
"\n",
45+
"# Deep heating parameters\n",
46+
"deep_P_dep = 1.0e6 # 10 bar deposition pressure\n",
47+
"deep_sigma_P = 1.0 # Width in log-pressure\n",
48+
"deep_efficiency = 0.02 # 2% of stellar flux\n",
49+
"\n",
50+
"# Expected additional flux\n",
51+
"F_deep_expected = deep_efficiency * instellation\n",
52+
"@printf(\"Expected deep heating flux: %.2f W/m²\\n\", F_deep_expected)"
53+
]
54+
},
55+
{
56+
"cell_type": "code",
57+
"execution_count": null,
58+
"id": "7cac6882",
59+
"metadata": {},
60+
"outputs": [],
61+
"source": [
62+
"# Function to create and run an atmosphere model\n",
63+
"function run_model(out_dir::String, deep_active::Bool)\n",
64+
" # Clean output directory\n",
65+
" rm(out_dir, force=true, recursive=true)\n",
66+
" mkpath(out_dir)\n",
67+
"\n",
68+
" # Create atmosphere\n",
69+
" atmos = AGNI.atmosphere.Atmos_t()\n",
70+
"\n",
71+
" # Composition: H2-dominated with some H2O\n",
72+
" mf_dict = Dict(\"H2\" => 0.85, \"He\" => 0.14, \"H2O\" => 0.01)\n",
73+
"\n",
74+
" # Setup atmosphere\n",
75+
" AGNI.atmosphere.setup!(atmos, ROOT_DIR, out_dir,\n",
76+
" spectral_file,\n",
77+
" instellation, 1.0, 0.0, 48.19, # instellation, s0_fact, albedo_b, zenith\n",
78+
" tmp_surf, gravity, radius,\n",
79+
" nlev, p_surf, p_top,\n",
80+
" mf_dict, \"\"; # mf_dict, mf_path\n",
81+
" flag_rayleigh=false,\n",
82+
" flag_gcontinuum=true,\n",
83+
" real_gas=false,\n",
84+
" thermo_functions=true)\n",
85+
"\n",
86+
" # Allocate\n",
87+
" AGNI.atmosphere.allocate!(atmos, joinpath(ROOT_DIR, star_file))\n",
88+
"\n",
89+
" # Set deep heating if requested\n",
90+
" if deep_active\n",
91+
" AGNI.atmosphere.set_deep_heating!(atmos,\n",
92+
" active=true,\n",
93+
" P_dep=deep_P_dep,\n",
94+
" sigma_P=deep_sigma_P,\n",
95+
" efficiency=deep_efficiency)\n",
96+
" end\n",
97+
"\n",
98+
" # Set initial temperature profile (isothermal + adiabatic)\n",
99+
" AGNI.setpt.isothermal!(atmos, tmp_surf)\n",
100+
"\n",
101+
" # Calculate layer properties\n",
102+
" AGNI.atmosphere.calc_layer_props!(atmos)\n",
103+
"\n",
104+
" # Calculate fluxes (no solve, just single evaluation)\n",
105+
" AGNI.energy.calc_fluxes!(atmos, true, false, true, false, true)\n",
106+
"\n",
107+
" return atmos\n",
108+
"end"
109+
]
110+
},
111+
{
112+
"cell_type": "code",
113+
"execution_count": null,
114+
"id": "cfb5b969",
115+
"metadata": {
116+
"vscode": {
117+
"languageId": "ruby"
118+
}
119+
},
120+
"outputs": [],
121+
"source": [
122+
"# Run model WITHOUT deep heating\n",
123+
"@info \"Running model without deep heating...\"\n",
124+
"out_dir_off = joinpath(OUT_DIR_BASE, \"no_deep_heating\")\n",
125+
"atmos_off = run_model(out_dir_off, false)\n",
126+
"@info \"Done\""
127+
]
128+
},
129+
{
130+
"cell_type": "code",
131+
"execution_count": null,
132+
"id": "6c6e6f80",
133+
"metadata": {},
134+
"outputs": [],
135+
"source": [
136+
"# Run model WITH deep heating\n",
137+
"@info \"Running model with deep heating...\"\n",
138+
"out_dir_on = joinpath(OUT_DIR_BASE, \"with_deep_heating\")\n",
139+
"atmos_on = run_model(out_dir_on, true)\n",
140+
"@info \"Done\""
141+
]
142+
},
143+
{
144+
"cell_type": "code",
145+
"execution_count": null,
146+
"id": "1ca11aba",
147+
"metadata": {},
148+
"outputs": [],
149+
"source": [
150+
"# Compare results\n",
151+
"println(\"=\"^60)\n",
152+
"println(\"COMPARISON OF RESULTS\")\n",
153+
"println(\"=\"^60)\n",
154+
"\n",
155+
"# TOA fluxes\n",
156+
"OLR_off = atmos_off.flux_u_lw[1]\n",
157+
"OLR_on = atmos_on.flux_u_lw[1]\n",
158+
"\n",
159+
"flux_tot_toa_off = atmos_off.flux_tot[1]\n",
160+
"flux_tot_toa_on = atmos_on.flux_tot[1]\n",
161+
"\n",
162+
"flux_tot_boa_off = atmos_off.flux_tot[end]\n",
163+
"flux_tot_boa_on = atmos_on.flux_tot[end]\n",
164+
"\n",
165+
"# Deep heating flux at BOA (should equal F_deep_expected)\n",
166+
"flux_deep_boa = atmos_on.flux_deep[end]\n",
167+
"\n",
168+
"@printf(\"\\nWithout Deep Heating:\\n\")\n",
169+
"@printf(\" OLR (LW up at TOA): %.2f W/m²\\n\", OLR_off)\n",
170+
"@printf(\" Total flux at TOA: %.2f W/m²\\n\", flux_tot_toa_off)\n",
171+
"@printf(\" Total flux at BOA: %.2f W/m²\\n\", flux_tot_boa_off)\n",
172+
"\n",
173+
"@printf(\"\\nWith Deep Heating (ε=%.2f):\\n\", deep_efficiency)\n",
174+
"@printf(\" OLR (LW up at TOA): %.2f W/m²\\n\", OLR_on)\n",
175+
"@printf(\" Total flux at TOA: %.2f W/m²\\n\", flux_tot_toa_on)\n",
176+
"@printf(\" Total flux at BOA: %.2f W/m²\\n\", flux_tot_boa_on)\n",
177+
"@printf(\" Deep heating flux (BOA): %.2f W/m²\\n\", flux_deep_boa)\n",
178+
"\n",
179+
"@printf(\"\\nDifferences:\\n\")\n",
180+
"@printf(\" ΔOLR: %.2f W/m²\\n\", OLR_on - OLR_off)\n",
181+
"@printf(\" ΔTotal flux (TOA): %.2f W/m²\\n\", flux_tot_toa_on - flux_tot_toa_off)\n",
182+
"@printf(\" ΔTotal flux (BOA): %.2f W/m²\\n\", flux_tot_boa_on - flux_tot_boa_off)\n",
183+
"\n",
184+
"@printf(\"\\nExpected deep heating flux: %.2f W/m²\\n\", F_deep_expected)\n",
185+
"@printf(\"Actual deep heating flux: %.2f W/m²\\n\", flux_deep_boa)\n",
186+
"@printf(\"Difference: %.4f W/m² (%.2f%%)\\n\",\n",
187+
" abs(flux_deep_boa - F_deep_expected),\n",
188+
" 100*abs(flux_deep_boa - F_deep_expected)/F_deep_expected)"
189+
]
190+
},
191+
{
192+
"cell_type": "code",
193+
"execution_count": null,
194+
"id": "f963a493",
195+
"metadata": {},
196+
"outputs": [],
197+
"source": [
198+
"# Plot temperature profiles\n",
199+
"p1 = plot(xlabel=\"Temperature [K]\", ylabel=\"Pressure [bar]\",\n",
200+
" yflip=true, yscale=:log10, legend=:topright,\n",
201+
" title=\"Temperature Profile\", size=(500, 400))\n",
202+
"\n",
203+
"plot!(p1, atmos_off.tmp, atmos_off.p .* 1e-5,\n",
204+
" label=\"No Deep Heating\", lw=2, color=:blue)\n",
205+
"plot!(p1, atmos_on.tmp, atmos_on.p .* 1e-5,\n",
206+
" label=\"With Deep Heating\", lw=2, color=:red, ls=:dash)\n",
207+
"\n",
208+
"display(p1)"
209+
]
210+
},
211+
{
212+
"cell_type": "code",
213+
"execution_count": null,
214+
"id": "1aa3bee6",
215+
"metadata": {},
216+
"outputs": [],
217+
"source": [
218+
"# Plot flux profiles\n",
219+
"p2 = plot(xlabel=\"Flux [W m⁻²]\", ylabel=\"Pressure [bar]\",\n",
220+
" yflip=true, yscale=:log10, legend=:topright,\n",
221+
" title=\"Total Flux Profile\", size=(500, 400))\n",
222+
"\n",
223+
"plot!(p2, atmos_off.flux_tot, atmos_off.pl .* 1e-5,\n",
224+
" label=\"Total (no deep)\", lw=2, color=:blue)\n",
225+
"plot!(p2, atmos_on.flux_tot, atmos_on.pl .* 1e-5,\n",
226+
" label=\"Total (with deep)\", lw=2, color=:red)\n",
227+
"plot!(p2, atmos_on.flux_deep, atmos_on.pl .* 1e-5,\n",
228+
" label=\"Deep Heating\", lw=2, color=:magenta, ls=:dash)\n",
229+
"\n",
230+
"# Mark deposition pressure\n",
231+
"hline!(p2, [deep_P_dep * 1e-5], ls=:dot, color=:gray, label=\"P_dep\")\n",
232+
"\n",
233+
"display(p2)"
234+
]
235+
},
236+
{
237+
"cell_type": "code",
238+
"execution_count": null,
239+
"id": "05050f1a",
240+
"metadata": {},
241+
"outputs": [],
242+
"source": [
243+
"# Plot deep heating flux gradient (dF/dP profile)\n",
244+
"# Calculate the flux gradient from the cumulative flux\n",
245+
"dF_deep = diff(atmos_on.flux_deep)\n",
246+
"dp = diff(atmos_on.pl)\n",
247+
"dF_dP = dF_deep ./ dp\n",
248+
"\n",
249+
"p3 = plot(xlabel=\"dF/dP [W m⁻² Pa⁻¹]\", ylabel=\"Pressure [bar]\",\n",
250+
" yflip=true, yscale=:log10, legend=:topright,\n",
251+
" title=\"Deep Heating Rate (dF/dP)\", size=(500, 400))\n",
252+
"\n",
253+
"plot!(p3, dF_dP, atmos_on.p .* 1e-5,\n",
254+
" label=\"dF_deep/dP\", lw=2, color=:magenta)\n",
255+
"\n",
256+
"# Mark deposition pressure\n",
257+
"hline!(p3, [deep_P_dep * 1e-5], ls=:dot, color=:gray, label=\"P_dep\")\n",
258+
"\n",
259+
"display(p3)"
260+
]
261+
},
262+
{
263+
"cell_type": "code",
264+
"execution_count": null,
265+
"id": "30699ab1",
266+
"metadata": {},
267+
"outputs": [],
268+
"source": [
269+
"# Combined plot with all flux components\n",
270+
"p4 = plot(xlabel=\"Flux [W m⁻²]\", ylabel=\"Pressure [bar]\",\n",
271+
" yflip=true, yscale=:log10, legend=:outerright,\n",
272+
" title=\"Flux Components (With Deep Heating)\", size=(700, 450))\n",
273+
"\n",
274+
"# Radiative fluxes\n",
275+
"plot!(p4, atmos_on.flux_u_lw, atmos_on.pl .* 1e-5,\n",
276+
" label=\"LW Up\", lw=2, color=:red)\n",
277+
"plot!(p4, -atmos_on.flux_d_lw, atmos_on.pl .* 1e-5,\n",
278+
" label=\"LW Down\", lw=2, color=:red, ls=:dash)\n",
279+
"plot!(p4, atmos_on.flux_u_sw, atmos_on.pl .* 1e-5,\n",
280+
" label=\"SW Up\", lw=2, color=:orange)\n",
281+
"plot!(p4, -atmos_on.flux_d_sw, atmos_on.pl .* 1e-5,\n",
282+
" label=\"SW Down\", lw=2, color=:orange, ls=:dash)\n",
283+
"\n",
284+
"# Convective flux\n",
285+
"plot!(p4, atmos_on.flux_cdry, atmos_on.pl .* 1e-5,\n",
286+
" label=\"Convection\", lw=2, color=:blue)\n",
287+
"\n",
288+
"# Deep heating flux\n",
289+
"plot!(p4, atmos_on.flux_deep, atmos_on.pl .* 1e-5,\n",
290+
" label=\"Deep Heating\", lw=3, color=:magenta)\n",
291+
"\n",
292+
"# Total flux\n",
293+
"plot!(p4, atmos_on.flux_tot, atmos_on.pl .* 1e-5,\n",
294+
" label=\"Total\", lw=2, color=:black)\n",
295+
"\n",
296+
"# Mark deposition pressure\n",
297+
"hline!(p4, [deep_P_dep * 1e-5], ls=:dot, color=:gray, lw=1, label=\"P_dep\")\n",
298+
"\n",
299+
"display(p4)"
300+
]
301+
},
302+
{
303+
"cell_type": "code",
304+
"execution_count": null,
305+
"id": "3cfe29c1",
306+
"metadata": {},
307+
"outputs": [],
308+
"source": [
309+
"# Save combined comparison plot\n",
310+
"p_combined = plot(p1, p2, p3, p4, layout=(2, 2), size=(1000, 800))\n",
311+
"savefig(p_combined, joinpath(OUT_DIR_BASE, \"deep_heating_comparison.png\"))\n",
312+
"@info \"Saved comparison plot to $(joinpath(OUT_DIR_BASE, \"deep_heating_comparison.png\"))\"\n",
313+
"display(p_combined)"
314+
]
315+
},
316+
{
317+
"cell_type": "code",
318+
"execution_count": null,
319+
"id": "62b1b230",
320+
"metadata": {},
321+
"outputs": [],
322+
"source": [
323+
"# Clean up\n",
324+
"AGNI.atmosphere.deallocate!(atmos_off)\n",
325+
"AGNI.atmosphere.deallocate!(atmos_on)\n",
326+
"@info \"Test complete!\""
327+
]
328+
}
329+
],
330+
"metadata": {
331+
"kernelspec": {
332+
"display_name": "venus",
333+
"language": "python",
334+
"name": "python3"
335+
},
336+
"language_info": {
337+
"name": "python",
338+
"version": "3.10.18"
339+
}
340+
},
341+
"nbformat": 4,
342+
"nbformat_minor": 5
343+
}

0 commit comments

Comments
 (0)