divertor.py

"""Module containing divertor routines"""

import math

import numpy as np

from process.core import constants
from process.core import process_output as po
from process.core.exceptions import ProcessValueError
from process.core.model import Model
from process.data_structure.divertor_variables import DivertorHeatLoadModel
from process.data_structure.physics_variables import DivertorNumberModels


class Divertor(Model):
    """Module containing divertor routines

    This module contains routines relevant for calculating the
    divertor parameters for a fusion power plant.
    """

    def __init__(self):
        self.outfile = constants.NOUT  # output file unit

    def output(self):
        """Output divertor information"""
        self.run(output=True)

    def run(self, output: bool = False):
        """Routine to call the divertor model

        This subroutine calls the divertor routine. This routine scales
        dimensions, powers and field levels which are used as input to
        the Harrison divertor model.

        Parameters
        ----------
        output :
            indicate whether output should be written to the output file, or not
        """
        self.data.divertor.deg_div_poloidal_plasma = self.single_divertor_angle
        self.data.fwbs.f_ster_div_single = (
            self.data.divertor.deg_div_poloidal_plasma / 360.0
        )

        self.data.fwbs.p_div_nuclear_heat_total_mw = self.incident_neutron_power(
            p_plasma_neutron_mw=self.data.physics.p_plasma_neutron_mw,
            f_ster_div_single=self.data.fwbs.f_ster_div_single,
            n_divertors=self.data.divertor.n_divertors,
        )

        self.data.fwbs.p_div_rad_total_mw = self.incident_radiation_power(
            p_plasma_rad_mw=self.data.physics.p_plasma_rad_mw,
            f_ster_div_single=self.data.fwbs.f_ster_div_single,
            n_divertors=self.data.divertor.n_divertors,
        )

        if (
            DivertorHeatLoadModel(self.data.divertor.i_div_heat_load)
            == DivertorHeatLoadModel.USER_INPUT
            and output
        ):
            po.ovarre(
                self.outfile,
                "Divertor heat load (MW/m²)",
                "(pflux_div_heat_load_mw)",
                self.data.divertor.pflux_div_heat_load_mw,
            )
            return
        if (
            DivertorHeatLoadModel(self.data.divertor.i_div_heat_load)
            == DivertorHeatLoadModel.PENG_CHAMBER
        ):
            self.divtart(
                self.data.physics.rmajor,
                self.data.physics.rminor,
                self.data.physics.triang,
                self.data.build.dr_fw_plasma_gap_inboard,
                self.data.build.dz_xpoint_divertor,
                self.data.physics.p_plasma_separatrix_mw,
                output=output,
                i_single_null=self.data.physics.i_single_null,
                dz_divertor=self.data.divertor.dz_divertor,
            )
            return
        if (
            DivertorHeatLoadModel(self.data.divertor.i_div_heat_load)
            == DivertorHeatLoadModel.WADE
        ):
            self.divwade(
                self.data.physics.rmajor,
                self.data.physics.rminor,
                self.data.physics.aspect,
                self.data.physics.b_plasma_toroidal_on_axis,
                self.data.physics.b_plasma_surface_poloidal_average,
                self.data.physics.p_plasma_separatrix_mw,
                self.data.divertor.f_div_flux_expansion,
                self.data.physics.nd_plasma_separatrix_electron,
                self.data.divertor.deg_div_field_plate,
                self.data.physics.rad_fraction_sol,
                self.data.physics.f_p_div_lower_separatrix,
                output=output,
            )
            return

    @property
    def single_divertor_angle(self):
        """
        Calculate the angle subtended by a single divertor.
        Angle is calculated as 180 degrees minus the inboard
        blanket poloidal angle, divided by 2 (for two divertors).
        """
        return (180.0 - self.data.blanket.deg_blkt_inboard_poloidal_plasma) / 2.0

    def divtart(
        self,
        rmajor: float,
        rminor: float,
        triang: float,
        dr_fw_plasma_gap_inboard: float,
        dz_xpoint_divertor: float,
        p_plasma_separatrix_mw: float,
        output: bool,
        i_single_null: int,
        dz_divertor: float,
    ) -> float:
        """Tight aspect ratio tokamak divertor model

        This method calculates the divertor heat load for a tight aspect
        ratio machine, assuming that the power is evenly distributed around the
        divertor chamber by the action of a gaseous target. Each divertor is
        modeled as approximately triangular in the R,Z plane.

        Parameters
        ----------
        rmajor : float
            Plasma major radius (m)
        rminor : float
            Plasma minor radius (m)
        triang : float
            Plasma triangularity
        dr_fw_plasma_gap_inboard : float
            Inboard scrape-off width (m)
        dz_xpoint_divertor : float
            Vertical distance from X-point to divertor (m)
        p_plasma_separatrix_mw : float
            Power to the divertor (MW)
        output : bool
            Indicates whether output should be written to the output file
        i_single_null : int
            1 for single null configuration, 0 for double null
        dz_divertor : float
            Vertical height of the divertor (m)

        Returns
        -------
        float
            Divertor heat load for a tight aspect ratio machine (MW/m2)

        Raises
        ------
        ProcessValueError
            If dz_xpoint_divertor is non-positive

        Notes
        -----
            - This model assumes a tight aspect ratio tokamak with a gaseous target
              divertor. The divertor chamber is modeled as triangular in the R,Z plane,
              and the heat load is calculated based on the total divertor surface area.
            - The method accounts for both single null and double null configurations.

        References
        ----------
            - Y.-K. M. Peng, J. B. Hicks, AEA Fusion, Culham (UK),
            "Engineering feasibility of tight aspect ratio Tokamak (spherical torus)
             reactors".
              1990. https://inis.iaea.org/records/ey2rf-dah04

            - Y.-K. M. Peng, J. B. Hicks,
            “Engineering feasibility of tight aspect ratio tokamak (spherical torus)
            reactors,”
              Osti.gov, 1991. https://www.osti.gov/biblio/1022679
              (accessed Mar. 24, 2025).
        """
        #  Thickness of centrepost + first wall at divertor height

        r1 = (
            rmajor
            - rminor * triang
            - 3.0e0 * dr_fw_plasma_gap_inboard
            + self.data.tfcoil.drtop
        )

        #  Outer radius of divertor region

        r2 = rmajor + rminor

        #  Angle of diagonal divertor plate from horizontal

        if dz_xpoint_divertor <= 0.0e0:
            raise ProcessValueError(
                "Non-positive dz_xpoint_divertor", dz_xpoint_divertor=dz_xpoint_divertor
            )

        theta = math.atan(dz_divertor / (r2 - r1))

        #  Vertical plate area

        a1 = 2.0e0 * np.pi * r1 * dz_divertor

        #  Horizontal plate area

        a2 = np.pi * (r2 * r2 - r1 * r1)

        #  Diagonal plate area

        a3 = a2 / (math.cos(theta) * math.cos(theta))

        #  Total divertor area

        # Single null case
        match DivertorNumberModels(i_single_null):
            case DivertorNumberModels.SINGLE_NULL:
                areadv = a1 + a2 + a3
            # Double null case
            case DivertorNumberModels.DOUBLE_NULL:
                areadv = 2.0 * (a1 + a2 + a3)

        if (
            DivertorHeatLoadModel(self.data.divertor.i_div_heat_load)
            == DivertorHeatLoadModel.PENG_CHAMBER
        ):
            self.data.divertor.pflux_div_heat_load_mw = p_plasma_separatrix_mw / areadv

        if (
            output
            and DivertorHeatLoadModel(self.data.divertor.i_div_heat_load)
            == DivertorHeatLoadModel.PENG_CHAMBER
        ):
            po.osubhd(self.outfile, "Divertor Heat Load")
            po.ocmmnt(self.outfile, "Assume an expanded divertor with a gaseous target")
            po.oblnkl(self.outfile)
            po.ovarre(
                self.outfile,
                "Power to the divertor (MW)",
                "(p_plasma_separatrix_mw.)",
                p_plasma_separatrix_mw,
            )
            po.ovarre(self.outfile, "Divertor surface area (m²)", "(areadv)", areadv)
            po.ovarre(
                self.outfile,
                "Divertor heat load (MW/m²)",
                "(pflux_div_heat_load_mw)",
                self.data.divertor.pflux_div_heat_load_mw,
            )

        elif output:
            po.osubhd(self.outfile, "Divertor Heat Load")
            po.ovarre(
                self.outfile,
                "Power to the divertor (MW)",
                "(p_plasma_separatrix_mw.)",
                p_plasma_separatrix_mw,
            )
            po.ovarre(
                self.outfile,
                "Divertor heat load (MW/m²)",
                "(pflux_div_heat_load_mw)",
                self.data.divertor.pflux_div_heat_load_mw,
            )
        return self.data.divertor.pflux_div_heat_load_mw

    def divwade(
        self,
        rmajor: float,
        rminor: float,
        aspect: float,
        b_plasma_toroidal_on_axis: float,
        b_plasma_poloidal_average: float,
        p_plasma_separatrix_mw: float,
        f_div_flux_expansion: float,
        nd_plasma_separatrix_electron: float,
        deg_div_field_plate: float,
        rad_fraction_sol: float,
        f_p_div_lower_separatrix: float,
        output: bool,
    ) -> float:
        """Divertor heat load model (Wade 2020)

        This subroutine calculates the divertor heat flux for any machine,
        with either a single null or double null configuration.
        It uses the Eich scaling (Eich et al. 2013) and spreading factor
        (Scarabosio et al. 2014)
        to calculate the SOL width.
        This is then used with a flux expansion factor to calculate
        the wetted area and then the heat flux.

        Parameters
        ----------
        rmajor : float
            plasma major radius (m)
        rminor : float
            plasma minor radius (m)
        aspect : float
            tokamak aspect ratio
        b_plasma_toroidal_on_axis : float
            toroidal field (T)
        b_plasma_poloidal_average : float
            poloidal field (T)
        p_plasma_separatrix_mw : float
            power to divertor (MW)
        f_div_flux_expansion : float
            plasma flux expansion in divertor
        nd_plasma_separatrix_electron : float
            electron density at separatrix (m-3)
        deg_div_field_plate : float
            field line angle wrt divertor target plate (degrees)
        rad_fraction_sol : float
            SOL radiation fraction
        f_p_div_lower_separatrix : float
            fraction of power to the lower divertor in double null configuration

        Returns
        -------
        float
            divertor heat load for a tight aspect ratio machine
        """
        # Radius on midplane
        r_omp = rmajor + rminor

        # B fields on midplane
        Bp_omp = -b_plasma_poloidal_average * rmajor / r_omp

        Bt_omp = -b_plasma_toroidal_on_axis * rmajor / r_omp

        # Eich scaling for lambda_q
        lambda_eich = (
            1.35
            * p_plasma_separatrix_mw**-0.02
            * rmajor**0.04
            * b_plasma_poloidal_average**-0.92
            * aspect**0.42
        )

        # Spreading factor
        spread_fact = (
            0.12
            * (nd_plasma_separatrix_electron / 1e19) ** -0.02
            * p_plasma_separatrix_mw**-0.21
            * rmajor**0.71
            * b_plasma_poloidal_average**-0.82
        )

        # SOL width
        lambda_int = lambda_eich + 1.64 * spread_fact

        # Flux angle on midplane
        alpha_mid = math.degrees(math.atan(Bp_omp / Bt_omp))

        # Flux angle in the divertor
        alpha_div = f_div_flux_expansion * alpha_mid

        # Tilt of the separatrix relative to the target in the poloidal plane
        theta_div = math.asin(
            (1 + 1 / alpha_div**2) * math.sin(math.radians(deg_div_field_plate))
        )

        # Wetted area
        area_wetted = (
            2 * np.pi * rmajor * lambda_int * f_div_flux_expansion * math.sin(theta_div)
        )

        # Divertor heat load
        hldiv_base = p_plasma_separatrix_mw * (1 - rad_fraction_sol) / area_wetted

        # For double null, calculate heat loads to upper and lower divertors
        # and use the highest
        if self.data.divertor.n_divertors == 2:
            hldiv_lower = f_p_div_lower_separatrix * hldiv_base
            hldiv_upper = (1.0 - f_p_div_lower_separatrix) * hldiv_base
            self.data.divertor.pflux_div_heat_load_mw = max(hldiv_lower, hldiv_upper)
        else:
            self.data.divertor.pflux_div_heat_load_mw = hldiv_base

        if output:
            po.osubhd(self.outfile, "Divertor Heat Load")
            po.ocmmnt(self.outfile, "Assume an expanded divertor with a gaseous target")
            po.oblnkl(self.outfile)
            po.ovarre(
                self.outfile,
                "Flux expansion",
                "(f_div_flux_expansion)",
                f_div_flux_expansion,
            )
            po.ovarre(
                self.outfile,
                "Field line angle wrt to target divertor plate (degrees)",
                "(deg_div_field_plate)",
                deg_div_field_plate,
            )
            po.ovarre(
                self.outfile,
                "Divertor heat load (MW/m²)",
                "(pflux_div_heat_load_mw)",
                self.data.divertor.pflux_div_heat_load_mw,
            )
        return self.data.divertor.pflux_div_heat_load_mw

    @staticmethod
    def incident_radiation_power(
        p_plasma_rad_mw: float,
        f_ster_div_single: float,
        n_divertors: int,
    ) -> float:
        """Calculates the total incident radiation power on the divertor box.

        Parameters
        ----------
        p_plasma_rad_mw : float
            Total plasma radiated power in megawatts (MW).
        f_ster_div_single : float
            Fraction of the solid angle subtended by a single divertor.
        n_divertors : int
            Number of divertors.

        Returns
        -------
        float
            Total incident radiation power on the divertor box in megawatts (MW).
        """
        return p_plasma_rad_mw * f_ster_div_single * n_divertors

    @staticmethod
    def incident_neutron_power(
        p_plasma_neutron_mw: float,
        f_ster_div_single: float,
        n_divertors: int,
    ) -> float:
        """Calculates the total incident neutron power on the divertor box.

        Parameters
        ----------
        p_plasma_neutron_mw : float
            Total plasma neutron power in megawatts (MW).
        f_ster_div_single : float
            Fraction of the solid angle subtended by a single divertor.
        n_divertors : int
            Number of divertors.

        Returns
        -------
        float
            Total incident radiation power on the divertor box in megawatts (MW).
        """
        return p_plasma_neutron_mw * f_ster_div_single * n_divertors


class LowerDivertor(Divertor):
    """Module containing lower divertor routines"""

    def run(self, output: bool):
        """Run the LowerDivertor routines"""
        super().run(output=output)

        self.data.divertor.p_div_lower_nuclear_heat_mw = self.incident_neutron_power(
            p_plasma_neutron_mw=self.data.physics.p_plasma_neutron_mw,
            f_ster_div_single=self.data.fwbs.f_ster_div_single,
            n_divertors=1,
        )

        self.data.divertor.p_div_lower_rad_mw = self.incident_radiation_power(
            p_plasma_rad_mw=self.data.physics.p_plasma_rad_mw,
            f_ster_div_single=self.data.fwbs.f_ster_div_single,
            n_divertors=1,
        )


class UpperDivertor(Divertor):
    """Module containing upper divertor routines"""

    def run(self, output: bool):
        """Run the UpperDivertor routine"""
        super().run(output=output)

        self.data.divertor.p_div_upper_nuclear_heat_mw = self.incident_neutron_power(
            p_plasma_neutron_mw=self.data.physics.p_plasma_neutron_mw,
            f_ster_div_single=self.data.fwbs.f_ster_div_single,
            n_divertors=1,
        )

        self.data.divertor.p_div_upper_rad_mw = self.incident_radiation_power(
            p_plasma_rad_mw=self.data.physics.p_plasma_rad_mw,
            f_ster_div_single=self.data.fwbs.f_ster_div_single,
            n_divertors=1,
        )