proteus.outgas.dummy¶

def calc_surface_pressures_dummy(dirs: dict, config: Config, hf_row: dict, *, initial: bool):
    """Compute volatile partitioning with parameterized model.

    Parameters
    ----------
    dirs : dict
        Directory paths.
    config : Config
        PROTEUS configuration.
    hf_row : dict
        Helpfile row (modified in place).
    initial : bool
        Whether this is an init-stage iteration. Only then is an empty O budget
        derived from the outgassed species; later an emptied budget stays empty.

    Notes
    -----
    The O in the atmosphere and melt is the stoichiometric O of the H2O, CO2 and
    SO2 built from H, C and S, and ``O_kg_total`` does not bound it. A user O
    budget below that O (``O_mode`` kg, ppmw or FeO_mantle_wt_pct), a budget a
    strip emptied, or an impact that delivers H, C or S without O leaves the
    atmosphere holding O that no budget holds: it sets ``P_surf`` and cannot
    escape, and an impact strips at most ``O_kg_total``.
    """
    Phi_global = float(hf_row['Phi_global'])
    gravity = float(hf_row['gravity'])
    R_int = float(hf_row['R_int'])
    if gravity <= 0 or R_int <= 0:
        raise ValueError(
            'Dummy outgassing needs a positive planet state: '
            f'gravity={gravity} m s-2, R_int={R_int} m'
        )
    area = 4.0 * np.pi * R_int**2

    # Save element totals before zeroing output keys (they are inputs)
    saved_element_kg = {}
    for element in element_list:
        saved_element_kg[element] = float(hf_row.get(f'{element}_kg_total', 0.0))

    # Initialize all expected keys to zero (M_atm is set by wrapper)
    for key in expected_keys():
        if key != 'M_atm':
            hf_row[key] = 0.0

    # Partition coefficient: fraction in atmosphere vs dissolved in melt.
    # Even at Phi=1 (fully molten), a fraction of volatiles is in the
    # atmosphere (solubility is finite). f_atm_floor ensures the
    # atmosphere is always non-empty so the coupling loop has a surface
    # pressure to work with.
    f_atm_floor = 0.1
    f_atm = f_atm_floor + (1.0 - f_atm_floor) * (1.0 - max(0.0, min(Phi_global, 1.0)))
    f_dissolved = 1.0 - f_atm

    # Convert elemental budgets to species masses
    species_kg_total = {}
    for element, (species, mass_frac) in _ELEMENT_TO_SPECIES.items():
        e_kg = saved_element_kg.get(element, 0.0)
        if e_kg > 0 and mass_frac > 0:
            species_kg_total[species] = e_kg / mass_frac  # kg of species

    # Partition each species between atmosphere and melt
    P_total = 0.0
    for species, kg_total in species_kg_total.items():
        kg_atm = f_atm * kg_total
        kg_liquid = f_dissolved * kg_total
        kg_solid = 0.0
        mmw = eval_gas_mmw(species)

        hf_row[f'{species}_kg_total'] = kg_total
        hf_row[f'{species}_kg_atm'] = kg_atm
        hf_row[f'{species}_kg_liquid'] = kg_liquid
        hf_row[f'{species}_kg_solid'] = kg_solid

        mol_total = kg_total / mmw
        hf_row[f'{species}_mol_total'] = mol_total
        hf_row[f'{species}_mol_atm'] = f_atm * mol_total
        hf_row[f'{species}_mol_liquid'] = f_dissolved * mol_total
        hf_row[f'{species}_mol_solid'] = 0.0

        # Partial pressure: P = m*g / A
        p_bar = kg_atm * gravity / area * 1e-5  # Pa -> bar
        hf_row[f'{species}_bar'] = p_bar
        P_total += p_bar

    hf_row['P_surf'] = P_total

    # VMRs from partial pressures
    for s in gas_list:
        if P_total > 0:
            hf_row[f'{s}_vmr'] = float(hf_row.get(f'{s}_bar', 0.0)) / P_total
        else:
            hf_row[f'{s}_vmr'] = 0.0

    # Mean molecular weight. An empty atmosphere (no volatile inventory)
    # has no meaningful MMW; store zero and let consumers check P_surf.
    if P_total > 0:
        mmw_sum = sum(hf_row.get(f'{s}_vmr', 0.0) * eval_gas_mmw(s) for s in gas_list)
        hf_row['atm_kg_per_mol'] = mmw_sum
    else:
        hf_row['atm_kg_per_mol'] = 0.0

    # Element reservoir masses (from species)
    _species_elements = {
        'H2O': {'H': 2 * 1.008 / 18.015, 'O': 15.999 / 18.015},
        'CO2': {'C': 12.011 / 44.009, 'O': 2 * 15.999 / 44.009},
        'N2': {'N': 1.0},
        'SO2': {'S': 32.060 / 64.058, 'O': 2 * 15.999 / 64.058},
    }
    for e in element_list:
        hf_row[f'{e}_kg_atm'] = 0.0
        hf_row[f'{e}_kg_liquid'] = 0.0
        hf_row[f'{e}_kg_solid'] = 0.0

    for species, kg_total in species_kg_total.items():
        elem_fracs = _species_elements.get(species, {})
        for e, frac in elem_fracs.items():
            hf_row[f'{e}_kg_atm'] += f_atm * kg_total * frac
            hf_row[f'{e}_kg_liquid'] += f_dissolved * kg_total * frac

    # Restore element totals (these are inputs, not outputs; the zeroing
    # above wiped them because expected_keys() includes _kg_total columns)
    for element, kg in saved_element_kg.items():
        if element != 'O':
            hf_row[f'{element}_kg_total'] = kg

    # Oxygen total: the saved budget is restored as given; an empty one is derived
    # from the outgassed species only in the init stage.
    saved_O = saved_element_kg.get('O', 0.0)
    if saved_O <= 0.0 and initial:
        saved_O = hf_row.get('O_kg_atm', 0.0) + hf_row.get('O_kg_liquid', 0.0)
    hf_row['O_kg_total'] = saved_O

    log.info(
        'Dummy outgas: P_surf=%.2f bar, Phi=%.3f, f_atm=%.3f, %d species',
        P_total,
        Phi_global,
        f_atm,
        len(species_kg_total),
    )