Skip to content

vulcan.agni

agni


Wrapper functions for coupling VULCAN to AGNI.

Imports
  • vulcan.config: Config class for handling configuration settings.
  • vulcan.paths: Paths to various files and directories used in VULCAN.
  • vulcan.phy_const: Physical constants used in VULCAN.
  • vulcan.store: AtmData and Variables classes for storing atmospheric data and variables.
  • vulcan.chem_funs: spec_list for the list of chemical species in VULCAN, gas_list for the list of gas species in AGNI.

_solve_energy(atmos, vulcan_cfg)

Use AGNI to solve for energy-conserving solution.

Source code in src/vulcan/agni.py
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
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
def _solve_energy(
    atmos,
    vulcan_cfg: Config,
):
    """Use AGNI to solve for energy-conserving solution."""

    # atmosphere solver plotting frequency
    modplot = 0
    if vulcan_cfg.use_live_plot:
        modplot = 1

    # tracking
    agni_success = False  # success?
    attempts = 0  # number of attempts so far

    # make attempts
    while not agni_success:
        attempts += 1
        log.info('Attempt %d' % attempts)

        # default parameters
        linesearch = 2
        easy_start = not bool(atmos.is_solved)
        dx_max = 300.0
        ls_increase = 1.04
        perturb_all = True
        max_steps = 70

        # try different solver parameters if struggling
        if attempts == 2:
            linesearch = 1
            dx_max *= 2.0
            ls_increase = 1.1

        log.debug('Solver parameters:')
        log.debug(
            '    ls_method=%d, easy_start=%s, dx_max=%.1f, ls_increase=%.2f'
            % (linesearch, str(easy_start), dx_max, ls_increase)
        )

        jl.AGNI.solver.ls_increase = float(ls_increase)

        # Try solving temperature profile
        agni_success = jl.AGNI.solver.solve_energy_b(
            atmos,
            sol_type=int(3),
            method=int(1),
            chem=False,
            conduct=False,
            convect=True,
            sens_heat=True,
            latent=False,
            rainout=True,
            max_steps=int(max_steps),
            max_runtime=900.0,
            conv_atol=float(vulcan_cfg.agni_atol),
            conv_rtol=float(vulcan_cfg.agni_rtol),
            ls_method=int(linesearch),
            dx_max=float(dx_max),
            easy_start=easy_start,
            perturb_all=perturb_all,
            save_frames=False,
            modplot=int(modplot),
            plot_jacobian=False,
        )

        # Model status check
        if agni_success:
            # success
            log.info('Attempt %d succeeded' % attempts)
            break
        else:
            # failure
            log.warning('Attempt %d failed' % attempts)

            # Max attempts
            if attempts >= 2:
                log.error('Maximum attempts when executing AGNI')
                break
    return atmos

_solve_once(atmos)

Use AGNI to solve radiative transfer with prescribed T(p) profile

Source code in src/vulcan/agni.py
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
def _solve_once(atmos):
    """Use AGNI to solve radiative transfer with prescribed T(p) profile"""

    #    dry convection
    jl.AGNI.setpt.dry_adiabat_b(atmos)

    #    temperature floor in stratosphere
    jl.AGNI.setpt.stratosphere_b(atmos, 0.5)

    # calculate convective flux only, to get Kzz profile
    jl.AGNI.energy.calc_fluxes_b(atmos, convective=True)

    # fill kzz values
    jl.AGNI.energy.fill_Kzz_b(atmos)

    return atmos

activate_julia(vulcan_cfg)

Source code in src/vulcan/agni.py
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
def activate_julia(vulcan_cfg: Config):

    # Check AGNI is in the right location
    if not os.path.isdir(paths.AGNI_DIR):
        raise FileNotFoundError(paths.AGNI_DIR)
    if not os.path.isfile(os.path.join(paths.AGNI_DIR, 'agni.jl')):
        raise FileNotFoundError(os.path.join(paths.AGNI_DIR, 'agni.jl'))

    # Activate environment for Julia
    log.info(f'Activating Julia environment {paths.AGNI_DIR}')
    jl.seval('using Pkg')
    jl.Pkg.activate(paths.AGNI_DIR)

    # Plotting configuration
    jl.seval('ENV["GKSwstype"] = "100"')
    jl.seval('using Plots')
    jl.seval('default(label=nothing, dpi=250)')

    # Import AGNI
    jl.seval('import AGNI')

    # Setup logging from AGNI
    #    This handle will be kept open throughout, so the file
    #    should not be deleted at runtime. However, it will be emptied when appropriate.
    verbosity = 2
    logpath = os.path.join(vulcan_cfg.output_dir, AGNI_LOGFILE_NAME)
    jl.AGNI.setup_logging(logpath, verbosity)

    log.debug("AGNI will log to '%s'" % logpath)

deallocate_atmos(atmos)

Deallocate atmosphere struct

Source code in src/vulcan/agni.py
171
172
173
174
175
def deallocate_atmos(atmos):
    """
    Deallocate atmosphere struct
    """
    jl.AGNI.atmosphere.deallocate_b(atmos)

init_agni_atmos(vulcan_cfg, atm, var)

Source code in src/vulcan/agni.py
 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
def init_agni_atmos(vulcan_cfg: Config, atm: AtmData, var: Variables):
    log.debug('New AGNI atmosphere')

    atmos = jl.AGNI.atmosphere.Atmos_t()
    succ = True

    # Stellar spectrum path (scaled to stellar surface)
    sflux_path = vulcan_cfg.sflux_file

    # Calculate instellation by scaling spectrum to TOA
    sflux_data = np.loadtxt(sflux_path).T  # erg/s/cm^2/nm
    sflux_integ = trapezoid(sflux_data[1], x=sflux_data[0])  # erg/s/cm^2
    sflux_integ *= (vulcan_cfg.r_star * r_sun) ** 2 / (vulcan_cfg.orbit_radius * au) ** 2
    sflux_integ *= 0.001  # W/m^2

    # Spectral file path
    try_spfile = os.path.join(vulcan_cfg.output_dir, 'runtime.sf')
    if os.path.exists(try_spfile):
        # exists => don't modify it
        input_sf = try_spfile
        input_star = ''
    elif vulcan_cfg.spectral_file == 'greygas':
        # grey gas solution
        input_sf = vulcan_cfg.spectral_file
        input_star = sflux_path
    else:
        # doesn't exist => AGNI will copy it + modify as required
        input_sf = os.path.join(paths.AGNI_DIR, vulcan_cfg.spectral_file)
        input_star = sflux_path

    # composition set initially well-mixed
    vol_dict = {}
    for g in gas_list:
        if '_' in g:
            continue
        vol_dict[g] = np.median(var.ymix[:, gas_list.index(g)])

    # set condensation
    condensates = []

    # Boundary pressures
    p_surf = atm.pico[0] * 1e-6  # convert to bar
    p_top = atm.pico[-1] * 1e-6  # convert to bar

    # FC working directory for AGNI
    fastchem_work = mkdtemp()

    # Setup struct
    succ = jl.AGNI.atmosphere.setup_b(
        atmos,
        paths.AGNI_DIR,
        vulcan_cfg.output_dir,
        input_sf,
        sflux_integ,
        vulcan_cfg.f_diurnal,
        0.0,
        vulcan_cfg.sl_angle * 180.0 / np.pi,  # convert radians to degrees
        vulcan_cfg.Tsurf_guess,
        vulcan_cfg.gs * 0.01,  # convert to SI, m/s^2
        vulcan_cfg.Rp * 0.01,  # convert to SI, m
        vulcan_cfg.agni_nlev,
        p_surf,
        p_top,
        vol_dict,
        '',
        flag_rayleigh=vulcan_cfg.use_rayleigh,
        flag_cloud=False,
        overlap_method='ee',
        albedo_s=vulcan_cfg.surf_albedo,
        surface_material='greybody',
        condensates=condensates,
        use_all_gases=False,
        fastchem_work=fastchem_work,
        real_gas=False,
        skin_d=0.01,
        skin_k=2.0,
        tmp_magma=vulcan_cfg.Tsurf_guess,
        tmp_floor=50.0,
    )
    if not succ:
        raise RuntimeError('Could not initialise AGNI atmos object')

    # Allocate arrays
    succ = jl.AGNI.atmosphere.allocate_b(atmos, input_star)

    if not succ:
        raise RuntimeError('Could not allocate AGNI atmos object')

    # Set temperature profile to initial guess
    tmp_top = 500.0
    tmp_top = min(tmp_top, vulcan_cfg.Tsurf_guess)
    log.debug('Initialised log-linear (top = %.2f K)' % tmp_top)
    jl.AGNI.setpt.loglinear_b(atmos, -0.5 * vulcan_cfg.Tsurf_guess)
    jl.AGNI.setpt.stratosphere_b(atmos, tmp_top)

    return atmos

run_agni(atmos, vulcan_cfg, atm, var)

Run AGNI atmosphere model.

Source code in src/vulcan/agni.py
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
def run_agni(atmos, vulcan_cfg: Config, atm: AtmData, var: Variables):
    """Run AGNI atmosphere model."""

    # Inform
    log.info('Running climate calculation...')

    # update gas compositions (interpolate from VULCAN to AGNI)
    for g in atmos.gas_names:
        p_vul = atm.pco[::-1] * 0.1  # convert to Pa
        g_vul = np.clip(var.ymix[::-1, gas_list.index(g)], 1e-30, 1.0)
        g_itp = PchipInterpolator(p_vul, np.log10(g_vul))

        atmos.gas_vmr[g][:] = 10 ** g_itp(atmos.p)[:]
        atmos.gas_ovmr[g][:] = atmos.gas_vmr[g][:]

    # jl.AGNI.plotting.plot_vmr(atmos,os.path.join(vulcan_cfg.plot_dir,"_agni_vmr.png"))

    # Run model
    if vulcan_cfg.solve_rce:
        log.info('Using nonlinear solver to conserve fluxes')
        atmos = _solve_energy(atmos, vulcan_cfg)
    else:
        log.info('Using prescribed temperature profile')
        atmos = _solve_once(atmos)

    # jl.AGNI.plotting.plot_pt(atmos,os.path.join(vulcan_cfg.plot_dir, "_agni_tp.png"))

    # Parse result (interpolate from AGNI to VULCAN)
    t_itp = PchipInterpolator(atmos.p, atmos.tmp)
    atm.Tco[:] = t_itp(atm.pco * 0.1)[:]

    log.info(spacer)
    return atmos