20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205 | def post_processing(config_params, id_mass=None, output_file=None, model_results=None):
"""Post-process model results by saving output data and plotting.
Parameters
----------
config_params : dict
Configuration parameters for the model.
id_mass : str or None
Identifier for the planet mass, used in output file naming.
output_file : str or None
Path to the output file for calculated mass and radius.
model_results : dict or None
Pre-computed solver output. When provided, ``post_processing``
skips its internal ``main()`` call (avoiding a redundant solve)
and writes output / plots directly from the supplied results.
Test helpers and benchmark drivers that already ran ``main()``
themselves should pass it here; the CLI keeps the implicit-call
path for back-compat.
"""
layer_eos_config = config_params['layer_eos_config']
solidus_id = config_params.get('rock_solidus', 'Stixrude14-solidus')
liquidus_id = config_params.get('rock_liquidus', 'Stixrude14-liquidus')
data_output_enabled = config_params['data_output_enabled']
plotting_enabled = config_params['plotting_enabled']
if model_results is None:
from .config import load_material_dictionaries
from .solver import main
model_results = main(
config_params,
material_dictionaries=load_material_dictionaries(),
melting_curves_functions=load_solidus_liquidus_functions(
layer_eos_config, solidus_id, liquidus_id
),
input_dir=os.path.join(get_zalmoxis_root(), 'input'),
)
# Extract results
radii = model_results['radii']
density = model_results['density']
gravity = model_results['gravity']
pressure = model_results['pressure']
temperature = model_results['temperature']
mass_enclosed = model_results['mass_enclosed']
cmb_mass = model_results['cmb_mass']
core_mantle_mass = model_results['core_mantle_mass']
total_time = model_results['total_time']
converged = model_results['converged']
converged_pressure = model_results['converged_pressure']
converged_density = model_results['converged_density']
converged_mass = model_results['converged_mass']
# Floor the CMB index at 1 so ``density[cmb_index - 1]`` (used for the
# core-side density log line) does not silently wrap to ``density[-1]``
# (the surface) when ``cmb_mass <= mass_enclosed[0]`` (e.g. coreless
# configs where np.argmax returns 0).
cmb_index = max(1, int(np.argmax(mass_enclosed >= cmb_mass)))
average_density = mass_enclosed[-1] / (4 / 3 * math.pi * radii[-1] ** 3)
# Check if mantle uses a Tdep EOS that needs external melting curves
# for phase detection. Unified PALEOS tables derive phases from the
# table itself and do not need (or support) get_Tdep_material().
mantle_str = layer_eos_config.get('mantle', '')
if mantle_str:
_mantle_mix = parse_layer_components(mantle_str)
uses_phase_detection = bool(set(_mantle_mix.components) & _NEEDS_MELTING_CURVES)
else:
uses_phase_detection = False
if uses_phase_detection:
mantle_pressures = pressure[cmb_index:]
mantle_temperatures = temperature[cmb_index:]
mantle_radii = radii[cmb_index:]
solidus_func, liquidus_func = load_solidus_liquidus_functions(
layer_eos_config, solidus_id, liquidus_id
)
mantle_phases = get_Tdep_material(
mantle_pressures, mantle_temperatures, solidus_func, liquidus_func
)
logger.info('Exoplanet Internal Structure Model Results:')
logger.info('----------------------------------------------------------------------')
logger.info(
f'Calculated Planet Mass: {mass_enclosed[-1]:.2e} kg or '
f'{mass_enclosed[-1] / earth_mass:.2f} Earth masses'
)
logger.info(
f'Calculated Planet Radius: {radii[-1]:.2e} m or '
f'{radii[-1] / earth_radius:.2f} Earth radii'
)
logger.info(f'Core Radius: {radii[cmb_index]:.2e} m')
logger.info(f'Mantle Density (at CMB): {density[cmb_index]:.2f} kg/m^3')
logger.info(f'Core Density (at CMB): {density[cmb_index - 1]:.2f} kg/m^3')
logger.info(f'Pressure at Core-Mantle Boundary (CMB): {pressure[cmb_index]:.2e} Pa')
logger.info(f'Pressure at Center: {pressure[0]:.2e} Pa')
logger.info(f'Average Density: {average_density:.2f} kg/m^3')
logger.info(f'CMB Mass Fraction: {mass_enclosed[cmb_index] / mass_enclosed[-1]:.3f}')
# ``core_mantle_mass = (cmf + mmf) * M_total`` is the cumulative mass at
# the top of the mantle layer, and ``mass_enclosed[cmb_index]`` is the
# core mass. The difference is therefore the mantle-layer mass alone,
# not ``core + mantle``. (For 2-layer configs with ``mmf = 0`` the
# mantle fills the remainder and this prints 0.000 by convention.)
logger.info(
f'Mantle Mass Fraction: '
f'{(core_mantle_mass - mass_enclosed[cmb_index]) / mass_enclosed[-1]:.3f}'
)
logger.info(f'Calculated Core Radius Fraction: {radii[cmb_index] / radii[-1]:.2f}')
logger.info(
f'Calculated Core+Mantle Radius Fraction: '
f'{(radii[np.argmax(mass_enclosed >= core_mantle_mass)] / radii[-1]):.2f}'
)
logger.info(f'Total Computation Time: {total_time:.2f} seconds')
logger.info(
f'Overall Convergence Status: {converged} with Pressure: {converged_pressure}, '
f'Density: {converged_density}, Mass: {converged_mass}'
)
if data_output_enabled:
output_data = np.column_stack(
(radii, density, gravity, pressure, temperature, mass_enclosed)
)
header = (
'Radius (m)\tDensity (kg/m^3)\tGravity (m/s^2)\t'
'Pressure (Pa)\tTemperature (K)\tMass Enclosed (kg)'
)
# Ensure the output directory exists. Local dev workstations
# typically have it from prior runs, but a fresh CI runner does
# not, and np.savetxt's writability probe fails there with
# FileNotFoundError. Creating it here keeps the writers
# independent of caller setup.
output_dir = os.path.join(get_zalmoxis_root(), 'output')
os.makedirs(output_dir, exist_ok=True)
if id_mass is None:
np.savetxt(
os.path.join(output_dir, 'planet_profile.txt'),
output_data,
header=header,
)
else:
np.savetxt(
os.path.join(output_dir, f'planet_profile{id_mass}.txt'),
output_data,
header=header,
)
if output_file is None:
output_file = os.path.join(
get_zalmoxis_root(), 'output', 'calculated_planet_mass_radius.txt'
)
if not os.path.exists(output_file):
header = 'Calculated Mass (kg)\tCalculated Radius (m)'
with open(output_file, 'w') as file:
file.write(header + '\n')
with open(output_file, 'a') as file:
file.write(f'{mass_enclosed[-1]}\t{radii[-1]}\n')
if plotting_enabled:
from .plots.plot_phase_vs_radius import plot_PT_with_phases
from .plots.plot_profiles import plot_planet_profile_single
plot_planet_profile_single(
radii,
density,
gravity,
pressure,
temperature,
radii[np.argmax(mass_enclosed >= cmb_mass)],
cmb_mass,
mass_enclosed[-1] / (4 / 3 * math.pi * radii[-1] ** 3),
mass_enclosed,
id_mass,
layer_eos_config=layer_eos_config,
)
if uses_phase_detection:
plot_PT_with_phases(
mantle_pressures,
mantle_temperatures,
mantle_radii,
mantle_phases,
radii[cmb_index],
)
|