diff --git a/CMakeLists.txt b/CMakeLists.txt
index c2c22a5153..f4e4744ccc 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -70,10 +70,6 @@ FILE(GLOB PROCESS_PYTHON_SRC_PATHS ${PROCESS_PYTHON_SRC_DIR}/*.py)
# Define interface source filenames to wrap
LIST(APPEND PROCESS_SRCS
numerics.f90
- constants.f90
- output.f90
- init_module.f90
- global_variables.f90
)
PREPROCESS()
diff --git a/process/availability.py b/process/availability.py
index ea1d65bff3..76c2a3161a 100644
--- a/process/availability.py
+++ b/process/availability.py
@@ -4,7 +4,7 @@
import numpy as np
from scipy.special import comb as combinations
-from process import fortran as ft
+from process import constants
from process import process_output as po
from process.data_structure import constraint_variables as ctv
from process.data_structure import cost_variables as cv
@@ -38,7 +38,7 @@ class Availability:
"""
def __init__(self) -> None:
- self.outfile = ft.constants.nout # output file unit
+ self.outfile = constants.NOUT # output file unit
def run(self, output: bool = False):
"""Run appropriate availability model
diff --git a/process/blanket_library.py b/process/blanket_library.py
index e9a8fb7296..e56429632b 100644
--- a/process/blanket_library.py
+++ b/process/blanket_library.py
@@ -7,6 +7,7 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.coolprop_interface import FluidProperties
from process.data_structure import (
@@ -19,7 +20,6 @@
primary_pumping_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import constants
logger = logging.getLogger(__name__)
@@ -40,7 +40,7 @@
class BlanketLibrary:
def __init__(self, fw) -> None:
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
self.fw = fw
diff --git a/process/build.py b/process/build.py
index 39cedf31b5..05fecf6e54 100644
--- a/process/build.py
+++ b/process/build.py
@@ -2,6 +2,7 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.blanket_library import dshellarea, eshellarea
from process.data_structure import (
@@ -15,15 +16,15 @@
tfcoil_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import constants, numerics
+from process.fortran import numerics
logger = logging.getLogger(__name__)
class Build:
def __init__(self):
- self.outfile = constants.nout
- self.mfile = constants.mfile
+ self.outfile = constants.NOUT
+ self.mfile = constants.MFILE
def run(self) -> None:
self.calculate_radial_build(output=False)
@@ -78,7 +79,7 @@ def calculate_beam_port_size(
# Have kept the single letter variable names to match the original code and documentation diagram.
radius_beam_tangency = f_radius_beam_tangency_rmajor * rmajor
- omega = constants.twopi / n_tf_coils
+ omega = constants.TWOPI / n_tf_coils
a = 0.5e0 * dx_tf_inboard_out_toroidal
try:
diff --git a/process/buildings.py b/process/buildings.py
index f0134124ed..f715f796c9 100644
--- a/process/buildings.py
+++ b/process/buildings.py
@@ -2,6 +2,7 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.data_structure import (
build_variables,
@@ -15,9 +16,6 @@
physics_variables,
tfcoil_variables,
)
-from process.fortran import (
- constants,
-)
logger = logging.getLogger(__name__)
@@ -34,7 +32,7 @@ def __init__(self) -> None:
This routine calls the buildings calculations.
"""
- self.outfile = constants.nout # output file unit
+ self.outfile = constants.NOUT # output file unit
def run(self, output: bool = False):
# Find TF coil radial positions
diff --git a/process/caller.py b/process/caller.py
index fb9d197163..282d95784c 100644
--- a/process/caller.py
+++ b/process/caller.py
@@ -14,7 +14,7 @@
from process.io.mfile import MFile
from process.iteration_variables import set_scaled_iteration_variable
from process.objectives import objective_function
-from process.utilities.f2py_string_patch import f2py_compatible_to_string
+from process.process_output import OutputFileManager
if TYPE_CHECKING:
from process.main import Models
@@ -132,16 +132,15 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
for _ in range(10):
# Divert OUT.DAT and MFILE.DAT output to scratch files for
# idempotence checking
- ft.init_module.open_idempotence_files()
+ OutputFileManager.open_idempotence_files()
self._call_models_once(xc)
# Write mfile
finalise(self.models, ifail)
# Extract data from intermediate idempotence-checking mfile
mfile_path = (
- f2py_compatible_to_string(ft.global_variables.output_prefix)
- + "IDEM_MFILE.DAT"
- )
+ data_structure.global_variables.output_prefix
+ ) + "IDEM_MFILE.DAT"
mfile = MFile(mfile_path)
# Create mfile dict of float values: only compare floats
mfile_data = {
@@ -176,7 +175,7 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
logger.debug("Mfiles idempotent, returning")
# Divert OUT.DAT and MFILE.DAT output back to original files
# now idempotence checking complete
- ft.init_module.close_idempotence_files()
+ OutputFileManager.close_idempotence_files()
# Write final output file and mfile
finalise(self.models, ifail)
return
@@ -202,7 +201,7 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
)
# Close idempotence files, write final output file and mfile
- ft.init_module.close_idempotence_files()
+ OutputFileManager.close_idempotence_files()
finalise(
self.models,
ifail,
@@ -213,7 +212,7 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int) -> None:
except Exception:
# If exception in model evaluations delete intermediate idempotence
# files to clean up
- ft.init_module.close_idempotence_files()
+ OutputFileManager.close_idempotence_files()
raise
def _call_models_once(self, xc: np.ndarray) -> None:
diff --git a/process/constants.py b/process/constants.py
new file mode 100644
index 0000000000..2c9c213651
--- /dev/null
+++ b/process/constants.py
@@ -0,0 +1,319 @@
+IOTTY = 6
+"""Standard output unit identifier"""
+
+NOUT = 11
+"""Output file unit identifier"""
+
+MFILE = 13
+"""Machine-optimised output file unit"""
+
+ELECTRON_CHARGE = 1.602176634e-19
+"""Electron / elementary charge [C]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?e|search_for=electron+charge
+"""
+
+ELECTRON_VOLT = ELECTRON_CHARGE
+"""Electron volt [J]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?evj|search_for=electron+volt
+"""
+
+KILOELECTRON_VOLT = ELECTRON_VOLT * 1e3
+"""Kiloelectron volt [kJ]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?evj|search_for=electron+volt
+"""
+
+ATOMIC_MASS_UNIT = 1.66053906892e-27
+"""Unified atomic mass unit [kg]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?u|search_for=atomic+mass+constant
+"""
+
+ELECTRON_MASS = 9.1093837139e-31
+"""Electron mass [kg]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?me|search_for=ELECTRON+MASS
+"""
+
+
+PROTON_MASS = 1.67262192595e-27
+"""Proton mass [kg]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?mp|search_for=PROTON+MASS
+"""
+
+M_PROTON_AMU = 1.0072764665789
+"""Proton mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?mpu
+"""
+
+M_PROTIUM_AMU = 1.00782503223
+"""Protium atomic mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=1
+"""
+
+DEUTERON_MASS = 3.3435837768e-27
+"""Deuteron mass [kg]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?md|search_for=DEUTERON+MASS
+"""
+
+M_DEUTERON_AMU = 2.013553212544
+"""Deuteron mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?mdu
+"""
+
+
+TRITON_MASS = 5.0073567512e-27
+"""Triton mass [kg]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?mt|search_for=TRITON+MASS
+"""
+
+M_TRITON_AMU = 3.01550071597
+"""Triton mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?art
+"""
+
+NEUTRON_MASS = 1.67492750056e-27
+"""Neutron mass [kg]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?mn|search_for=NEUTRON+MASS
+"""
+
+ALPHA_MASS = 6.6446573450e-27
+"""Alpha particle mass [kg]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?mal|search_for=alpha+mass
+"""
+
+M_ALPHA_AMU = 4.001506179129
+"""Alpha particle mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?malu
+"""
+
+M_HELIUM_AMU = 4.002602
+"""Average Helium atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=He
+"""
+
+HELION_MASS = 5.0064127862e-27
+"""Helion (3He) mass [kg]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?mh|search_for=HELION
+"""
+
+M_HELION_AMU = 3.014932246932
+"""Helion (3He) mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?arh
+"""
+
+M_BERYLLIUM_AMU = 9.0121831
+"""Beryllium atom (4Be) mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Be
+"""
+
+M_CARBON_AMU = 12.0096
+"""Average Carbon atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=C
+"""
+
+M_NITROGEN_AMU = 14.00643
+"""Average Nitrogen atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=N
+"""
+
+M_OXYGEN_AMU = 15.99903
+"""Average Oxygen atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=O
+"""
+
+M_NEON_AMU = 20.1797
+"""Average Neon atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Ne
+"""
+
+M_SILICON_AMU = 28.084
+"""Average Silicon atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Si
+"""
+
+M_ARGON_AMU = 39.948
+"""Average Argon atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Ar
+"""
+
+M_IRON_AMU = 55.845
+"""Average Iron atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Fe
+"""
+
+M_NICKEL_AMU = 58.6934
+"""Average Nickel atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Ni
+"""
+
+M_KRYPTON_AMU = 83.798
+"""Average Krypton atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Kr
+"""
+
+M_XENON_AMU = 131.293
+"""Average Xenon atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Xe
+"""
+
+M_TUNGSTEN_AMU = 183.84
+"""Average Tungsten atom mass [amu]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=W
+"""
+
+SPEED_LIGHT = 299792458
+"""Speed of light in vacuum (c) [m/s]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?c|search_for=light
+"""
+
+D_T_ENERGY = (
+ (DEUTERON_MASS + TRITON_MASS) - (ALPHA_MASS + NEUTRON_MASS)
+) * SPEED_LIGHT**2
+"""Deuterium - Tritium reaction energy [J]
+Find the mass difference in the reactancts and products of the D-T reaction
+Multiply by the speed of light squared to get the energy released
+"""
+
+D_HELIUM_ENERGY = (
+ (DEUTERON_MASS + HELION_MASS) - (ALPHA_MASS + PROTON_MASS)
+) * SPEED_LIGHT**2
+"""Deuterium - Helion (3He) reaction energy [J]
+Find the mass difference in the reactancts and products of the D-3He reaction
+Multiply by the speed of light squared to get the energy released
+"""
+
+DD_HELIUM_ENERGY = (
+ (DEUTERON_MASS + DEUTERON_MASS) - (HELION_MASS + NEUTRON_MASS)
+) * SPEED_LIGHT**2
+"""Deuterium - Deuterium (3He producing) reaction energy [J]
+Find the mass difference in the reactancts and products of the D-D reaction
+Multiply by the speed of light squared to get the energy released
+"""
+
+DD_TRITON_ENERGY = (
+ (DEUTERON_MASS + DEUTERON_MASS) - (TRITON_MASS + PROTON_MASS)
+) * SPEED_LIGHT**2
+"""Deuterium - Deuterium (Triton producing) reaction energy [J]
+Find the mass difference in the reactancts and products of the D-D reaction
+Multiply by the speed of light squared to get the energy released
+"""
+
+DT_NEUTRON_ENERGY_FRACTION = ALPHA_MASS / (NEUTRON_MASS + ALPHA_MASS)
+"""Deuterium - Tritium reaction energy fraction carried by neutron
+Assuming centre of mass frame as the momenta of the fusion products exceed
+those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
+Roughly 79.867% of the energy is carried by the neutron
+"""
+
+DT_ALPHA_ENERGY = (1.0 - DT_NEUTRON_ENERGY_FRACTION) * D_T_ENERGY
+"""Deuterium - Tritium reaction energy carried by alpha particle neutron [J]
+Assuming centre of mass frame as the momenta of the fusion products exceed
+those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
+Roughly 3.5 MeV of the energy is carried by the alpha particle
+"""
+
+DD_NEUTRON_ENERGY_FRACTION = HELION_MASS / (NEUTRON_MASS + HELION_MASS)
+"""Deuterium - Deuterium (3He producing) reaction energy fraction carried by neutron
+Assuming centre of mass frame as the momenta of the fusion products exceed
+those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
+Roughly 74.935% of the energy is carried by the neutron
+"""
+
+DD_PROTON_ENERGY_FRACTION = TRITON_MASS / (PROTON_MASS + TRITON_MASS)
+"""Deuterium - Deuterium (Triton producing) reaction energy fraction carried by proton
+Assuming centre of mass frame as the momenta of the fusion products exceed
+those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
+Roughly 74.960% of the energy is carried by the proton
+"""
+
+DHELIUM_PROTON_ENERGY_FRACTION = ALPHA_MASS / (PROTON_MASS + ALPHA_MASS)
+"""Deuterium - Helion (3He) reaction energy fraction carried by proton
+Assuming centre of mass frame as the momenta of the fusion products exceed
+those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
+Roughly 79.889% of the energy is carried by the proton
+"""
+
+DEN_TUNGSTEN = 19250.0
+"""Density of Tungsten [kg/m3]"""
+
+TEMP_ROOM = 293.15
+""" Room temperature in Kelvin
+Assume the room is at 20 degrees Celsius
+"""
+
+PI = 3.1415926535897932
+
+RMU0 = 1.256637062e-6
+"""permeability of free space [H/m]"""
+
+TWOPI = 6.2831853071795862
+
+UMASS = 1.660538921e-27
+"""unified atomic mass unit [kg]"""
+
+EPSILON0 = 8.85418781e-12
+"""permittivity of free space [Farad/m]"""
+
+CPH2O = 4180.0
+"""specific heat capacity of water (J/kg/K)"""
+
+den_copper: float = None
+"""density of copper (kg/m3)"""
+
+den_aluminium: float = None
+"""density of aluminium (kg/m3)"""
+
+DENH2O = 985.0
+"""density of water (kg/m3)"""
+
+K_COPPER = 330.0
+"""Copper thermal conductivity (W/m/K)"""
+
+KH2O = 0.651
+"""thermal conductivity of water (W/m/K)"""
+
+MUH2O = 4.71e-4
+"""water dynamic viscosity (kg/m/s)"""
+
+N_DAY_YEAR = 365.2425
+"""Average number of days in a year"""
+
+ACCELERATION_GRAVITY = 9.81
+"""Acceleration due to gravity [m/s2]"""
+
+
+def init_constants():
+ global den_copper
+ global den_aluminium
+
+ den_copper = 8900.0
+ den_aluminium = 2700.0
diff --git a/process/constraints.py b/process/constraints.py
index c91357b34e..3b32e3d85e 100644
--- a/process/constraints.py
+++ b/process/constraints.py
@@ -7,6 +7,7 @@
import process.data_structure as data_structure
import process.fortran as fortran
+from process import constants
from process.exceptions import ProcessError, ProcessValueError
ConstraintSymbolType = Literal["=", ">=", "<="]
@@ -142,8 +143,8 @@ def constraint_equation_1():
data_structure.physics_variables.beta_fast_alpha
+ data_structure.physics_variables.beta_beam
+ 2.0e3
- * fortran.constants.rmu0
- * fortran.constants.electron_charge
+ * constants.RMU0
+ * constants.ELECTRON_CHARGE
* (
data_structure.physics_variables.dene
* data_structure.physics_variables.ten
diff --git a/process/costs.py b/process/costs.py
index f3718d2cf9..07f5e41258 100644
--- a/process/costs.py
+++ b/process/costs.py
@@ -1,5 +1,6 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.data_structure import (
build_variables,
@@ -20,14 +21,11 @@
vacuum_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import (
- constants,
-)
class Costs:
def __init__(self):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def run(self):
"""
@@ -1692,7 +1690,7 @@ def acc2222(self):
for i in range(pfcoil_variables.n_cs_pf_coils):
pfwndl = (
pfwndl
- + constants.twopi
+ + constants.TWOPI
* pfcoil_variables.r_pf_coil_middle[i]
* pfcoil_variables.n_pf_coil_turns[i]
)
@@ -1789,7 +1787,7 @@ def acc2222(self):
cost_variables.c22221 = cost_variables.c22221 + (
1.0e-6
- * constants.twopi
+ * constants.TWOPI
* pfcoil_variables.r_pf_coil_middle[i]
* pfcoil_variables.n_pf_coil_turns[i]
* cpfconpm
@@ -1867,7 +1865,7 @@ def acc2222(self):
cost_variables.c22221 = cost_variables.c22221 + (
1.0e-6
- * constants.twopi
+ * constants.TWOPI
* pfcoil_variables.r_pf_coil_middle[pfcoil_variables.n_cs_pf_coils - 1]
* pfcoil_variables.n_pf_coil_turns[pfcoil_variables.n_cs_pf_coils - 1]
* cpfconpm
@@ -2459,7 +2457,7 @@ def acc2272(self):
2.0e0
* physics_variables.rndfuel
* physics_variables.m_fuel_amu
- * constants.umass
+ * constants.UMASS
* 1000.0e0
* 86400.0e0
)
@@ -2470,7 +2468,7 @@ def acc2272(self):
* 3.0e0
* 1.67e-27
* 1.0e3
- / (constants.electron_volt * 17.6e6 * ife_variables.fburn)
+ / (constants.ELECTRON_VOLT * 17.6e6 * ife_variables.fburn)
)
physics_variables.wtgpd = targtm * ife_variables.reprat * 86400.0e0
@@ -2860,14 +2858,14 @@ def coelc(self):
kwhpy = (
1.0e3
* heat_transport_variables.p_plant_electric_net_mw
- * (24.0e0 * constants.n_day_year)
+ * (24.0e0 * constants.N_DAY_YEAR)
* cost_variables.cfactr
)
else:
kwhpy = (
1.0e3
* heat_transport_variables.p_plant_electric_net_mw
- * (24.0e0 * constants.n_day_year)
+ * (24.0e0 * constants.N_DAY_YEAR)
* cost_variables.cfactr
* times_variables.t_burn
/ times_variables.t_cycle
@@ -3079,7 +3077,7 @@ def coelc(self):
* physics_variables.wtgpd
* 1.0e-3
* cost_variables.uche3
- * constants.n_day_year
+ * constants.N_DAY_YEAR
* cost_variables.cfactr
)
else:
diff --git a/process/costs_2015.py b/process/costs_2015.py
index bd67819b67..6f2a6c1e3f 100644
--- a/process/costs_2015.py
+++ b/process/costs_2015.py
@@ -2,6 +2,7 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.data_structure import (
build_variables,
@@ -9,23 +10,20 @@
cost_variables,
current_drive_variables,
fwbs_variables,
+ global_variables,
heat_transport_variables,
pf_power_variables,
pfcoil_variables,
physics_variables,
tfcoil_variables,
)
-from process.fortran import (
- constants,
- global_variables,
-)
logger = logging.getLogger(__name__)
class Costs2015:
def __init__(self):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def run(self):
"""
@@ -296,7 +294,7 @@ def calc_fwbs_costs(self):
cost_2015_variables.s_k[24] = (
build_variables.a_fw_total
* fwbs_variables.fw_armour_thickness
- * constants.den_tungsten
+ * constants.DEN_TUNGSTEN
)
cost_2015_variables.s_kref[24] = 29000.0e0
cost_2015_variables.s_cost[24] = (
diff --git a/process/cryostat.py b/process/cryostat.py
index 66d4aec157..ab4951d82a 100644
--- a/process/cryostat.py
+++ b/process/cryostat.py
@@ -1,5 +1,6 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.data_structure import (
blanket_library,
@@ -8,14 +9,11 @@
fwbs_variables,
pfcoil_variables,
)
-from process.fortran import (
- constants,
-)
class Cryostat:
def __init__(self) -> None:
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def run(self) -> None:
"""Run the cryostat calculations.
diff --git a/process/cs_fatigue.py b/process/cs_fatigue.py
index 5738803a47..d7840f1055 100644
--- a/process/cs_fatigue.py
+++ b/process/cs_fatigue.py
@@ -2,12 +2,12 @@
from numba import njit
import process.data_structure as data_structure
-from process.fortran import constants
+from process import constants
class CsFatigue:
def __init__(self):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def ncycle(
self,
diff --git a/process/current_drive.py b/process/current_drive.py
index 713df0e4bd..4ad1707ae5 100644
--- a/process/current_drive.py
+++ b/process/current_drive.py
@@ -2,6 +2,7 @@
import numpy as np
+from process import constants
from process import (
process_output as po,
)
@@ -11,7 +12,6 @@
physics_variables,
)
from process.exceptions import ProcessError, ProcessValueError
-from process.fortran import constants
from process.plasma_profiles import PlasmaProfile
logger = logging.getLogger(__name__)
@@ -19,7 +19,7 @@
class NeutralBeam:
def __init__(self, plasma_profile: PlasmaProfile):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
self.plasma_profile = plasma_profile
def iternb(self):
@@ -481,7 +481,7 @@ def cfnbi(self, afast, efast, te, ne, _nd, _nt, zeffai, xlmbda):
atmdt = 2.5
# atmt = 3.0
c = 3.0e8
- me = constants.electron_mass
+ me = constants.ELECTRON_MASS
# zd = 1.0
# zt = 1.0
@@ -494,17 +494,17 @@ def cfnbi(self, afast, efast, te, ne, _nd, _nt, zeffai, xlmbda):
xlmbdai = self.xlmbdabi(afast, atmdt, efast, te, ne)
sumln = zeffai * xlmbdai / xlmbda
xlnrat = (
- 3.0e0 * np.sqrt(np.pi) / 4.0e0 * me / constants.proton_mass * sumln
+ 3.0e0 * np.sqrt(np.pi) / 4.0e0 * me / constants.PROTON_MASS * sumln
) ** (2.0e0 / 3.0e0)
ve = c * np.sqrt(2.0e0 * te / 511.0e0)
ecritfi = (
afast
- * constants.proton_mass
+ * constants.PROTON_MASS
* ve
* ve
* xlnrat
- / (2.0e0 * constants.electron_charge * 1.0e3)
+ / (2.0e0 * constants.ELECTRON_CHARGE * 1.0e3)
)
x = np.sqrt(efast / ecritfi)
@@ -536,7 +536,7 @@ def xlmbdabi(self, mb, mth, eb, t, nelec):
class ElectronCyclotron:
def __init__(self, plasma_profile: PlasmaProfile):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
self.plasma_profile = plasma_profile
def culecd(self):
@@ -626,7 +626,7 @@ def eccdef(self, tlocal, epsloc, zlocal, cosang, coulog):
ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990
"""
mcsq = (
- constants.electron_mass * 2.9979e8**2 / (1.0e3 * constants.electron_volt)
+ constants.ELECTRON_MASS * 2.9979e8**2 / (1.0e3 * constants.ELECTRON_VOLT)
) # keV
f = 16.0e0 * (tlocal / mcsq) ** 2
@@ -748,7 +748,7 @@ def electron_cyclotron_freethy(
"""
# Cyclotron frequency
- fc = 1 / (2 * np.pi) * constants.electron_charge * bt / constants.electron_mass
+ fc = 1 / (2 * np.pi) * constants.ELECTRON_CHARGE * bt / constants.ELECTRON_MASS
# Plasma frequency
fp = (
@@ -756,8 +756,8 @@ def electron_cyclotron_freethy(
/ (2 * np.pi)
* np.sqrt(
(dene / 1.0e19)
- * constants.electron_charge**2
- / (constants.electron_mass * constants.epsilon0)
+ * constants.ELECTRON_CHARGE**2
+ / (constants.ELECTRON_MASS * constants.EPSILON0)
)
)
@@ -848,7 +848,7 @@ def legend(self, zlocal, arg):
class IonCyclotron:
def __init__(self, plasma_profile: PlasmaProfile):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
self.plasma_profile = plasma_profile
def ion_cyclotron_ipdg89(
@@ -889,7 +889,7 @@ def ion_cyclotron_ipdg89(
class ElectronBernstein:
def __init__(self, plasma_profile: PlasmaProfile):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
self.plasma_profile = plasma_profile
def electron_bernstein_freethy(
@@ -946,9 +946,9 @@ def electron_bernstein_freethy(
1.0e0
/ (2.0e0 * np.pi)
* n_ecrh_harmonic
- * constants.electron_charge
+ * constants.ELECTRON_CHARGE
* bt
- / constants.electron_mass
+ / constants.ELECTRON_MASS
)
fp = (
@@ -956,8 +956,8 @@ def electron_bernstein_freethy(
/ (2.0e0 * np.pi)
* np.sqrt(
dene20
- * constants.electron_charge**2
- / (constants.electron_mass * constants.epsilon0)
+ * constants.ELECTRON_CHARGE**2
+ / (constants.ELECTRON_MASS * constants.EPSILON0)
)
)
@@ -968,7 +968,7 @@ def electron_bernstein_freethy(
class LowerHybrid:
def __init__(self, plasma_profile: PlasmaProfile):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
self.plasma_profile = plasma_profile
def cullhy(self):
@@ -1241,7 +1241,7 @@ def __init__(
neutral_beam: NeutralBeam,
electron_bernstein: ElectronBernstein,
):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
self.plasma_profile = plasma_profile
self.electron_cyclotron = electron_cyclotron
self.ion_cyclotron = ion_cyclotron
diff --git a/process/data_structure/__init__.py b/process/data_structure/__init__.py
index 9f95874990..6ce2b43649 100644
--- a/process/data_structure/__init__.py
+++ b/process/data_structure/__init__.py
@@ -11,6 +11,7 @@
dcll_variables,
divertor_variables,
fwbs_variables,
+ global_variables,
heat_transport_variables,
ife_variables,
impurity_radiation_module,
@@ -47,6 +48,7 @@
"dcll_variables",
"divertor_variables",
"fwbs_variables",
+ "global_variables",
"heat_transport_variables",
"ife_variables",
"impurity_radiation_module",
diff --git a/process/data_structure/global_variables.py b/process/data_structure/global_variables.py
new file mode 100644
index 0000000000..e961ce3a66
--- /dev/null
+++ b/process/data_structure/global_variables.py
@@ -0,0 +1,65 @@
+icase: str = None
+"""Power plant type"""
+
+runtitle: str
+"""A short descriptive title for the run"""
+
+run_tests: int = None
+
+verbose: int = None
+
+maxcal: int = None
+"""Maximum number of solver iterations"""
+
+fileprefix: str = None
+"""Input file path prefix"""
+
+output_prefix: str = None
+"""Output file path prefix"""
+
+xlabel: str = None
+"""Scan parameter description label"""
+
+vlabel: str = None
+"""Scan value name label"""
+
+xlabel_2: str = None
+"""Scan parameter description label (2nd dimension)"""
+
+vlabel_2: str = None
+"""Scan value name label (2nd dimension)"""
+
+iscan_global: int = None
+
+convergence_parameter: int = None
+"""VMCON convergence parameter 'sum'"""
+
+
+def init_global_variables():
+ global icase
+ global runtitle
+ global run_tests
+ global verbose
+ global maxcal
+ global fileprefix
+ global output_prefix
+ global xlabel
+ global vlabel
+ global xlabel_2
+ global vlabel_2
+ global iscan_global
+ global convergence_parameter
+
+ icase = "Steady-state tokamak model"
+ runtitle = "Run Title (change this line using input variable 'runtitle')"
+ run_tests = 0
+ verbose = 0
+ maxcal = 200
+ fileprefix = ""
+ output_prefix = ""
+ xlabel = ""
+ vlabel = ""
+ xlabel_2 = ""
+ vlabel_2 = ""
+ iscan_global = 0
+ convergence_parameter = 0.0
diff --git a/process/dcll.py b/process/dcll.py
index e02831ec23..636bb90bc8 100644
--- a/process/dcll.py
+++ b/process/dcll.py
@@ -1,4 +1,5 @@
import process.blanket_library as blanket_library
+from process import constants
from process import (
process_output as po,
)
@@ -12,9 +13,6 @@
physics_variables,
primary_pumping_variables,
)
-from process.fortran import (
- constants,
-)
class DCLL(BlanketLibrary):
@@ -699,7 +697,7 @@ def dcll_masses(self, output: bool):
)
# First wall armour mass (kg)
fwbs_variables.fw_armour_mass = (
- fwbs_variables.fw_armour_vol * constants.den_tungsten
+ fwbs_variables.fw_armour_vol * constants.DEN_TUNGSTEN
)
# Total mass of blanket
diff --git a/process/divertor.py b/process/divertor.py
index 71a906be52..7713c873bc 100644
--- a/process/divertor.py
+++ b/process/divertor.py
@@ -1,12 +1,12 @@
import math
+from process import constants
from process import process_output as po
from process.data_structure import build_variables as bv
from process.data_structure import divertor_variables as dv
from process.data_structure import physics_variables as pv
from process.data_structure import tfcoil_variables as tfv
from process.exceptions import ProcessValueError
-from process.fortran import constants
class Divertor:
@@ -18,7 +18,7 @@ class Divertor:
"""
def __init__(self) -> None:
- self.outfile = constants.nout # output file unit
+ self.outfile = constants.NOUT # output file unit
def run(self, output: bool) -> None:
"""Routine to call the divertor model
@@ -153,11 +153,11 @@ def divtart(
# Vertical plate area
- a1 = 2.0e0 * constants.pi * r1 * dz_divertor
+ a1 = 2.0e0 * constants.PI * r1 * dz_divertor
# Horizontal plate area
- a2 = constants.pi * (r2 * r2 - r1 * r1)
+ a2 = constants.PI * (r2 * r2 - r1 * r1)
# Diagonal plate area
@@ -324,7 +324,7 @@ def divwade(
# Wetted area
area_wetted = (
2
- * constants.pi
+ * constants.PI
* rmajor
* lambda_int
* f_div_flux_expansion
diff --git a/process/evaluators.py b/process/evaluators.py
index 9eb5da7779..ab65e14379 100644
--- a/process/evaluators.py
+++ b/process/evaluators.py
@@ -5,9 +5,9 @@
from process.caller import Caller
from process.data_structure import cost_variables as cv
+from process.data_structure import global_variables as gv
from process.data_structure import physics_variables as pv
from process.data_structure import times_variables as tv
-from process.fortran import global_variables as gv
from process.fortran import numerics
logger = logging.getLogger(__name__)
diff --git a/process/final.py b/process/final.py
index d0d4e495a2..062076ae11 100644
--- a/process/final.py
+++ b/process/final.py
@@ -3,14 +3,10 @@
from tabulate import tabulate
import process.constraints as constraints
+from process import constants
from process import output as op
-from process import (
- process_output as po,
-)
-from process.fortran import (
- constants,
- numerics,
-)
+from process import process_output as po
+from process.fortran import numerics
from process.objectives import objective_function
from process.utilities.f2py_string_patch import f2py_compatible_to_string
@@ -28,9 +24,9 @@ def finalise(models, ifail: int, non_idempotent_msg: None | str = None):
:type non_idempotent_msg: None | str, optional
"""
if ifail == 1:
- po.oheadr(constants.nout, "Final Feasible Point")
+ po.oheadr(constants.NOUT, "Final Feasible Point")
else:
- po.oheadr(constants.nout, "Final UNFEASIBLE Point")
+ po.oheadr(constants.NOUT, "Final UNFEASIBLE Point")
# Output relevant to no optimisation
if numerics.ioptimz == -2:
@@ -38,23 +34,23 @@ def finalise(models, ifail: int, non_idempotent_msg: None | str = None):
# Print non-idempotence warning to OUT.DAT only
if non_idempotent_msg:
- po.oheadr(constants.nout, "NON-IDEMPOTENT VARIABLES")
- po.ocmmnt(constants.nout, non_idempotent_msg)
+ po.oheadr(constants.NOUT, "NON-IDEMPOTENT VARIABLES")
+ po.ocmmnt(constants.NOUT, non_idempotent_msg)
# Write output to OUT.DAT and MFILE.DAT
- op.write(models, constants.nout)
+ op.write(models, constants.NOUT)
def output_evaluation():
"""Write output for an evaluation run of PROCESS"""
- po.oheadr(constants.nout, "Numerics")
- po.ocmmnt(constants.nout, "PROCESS has performed an evaluation run.")
- po.oblnkl(constants.nout)
+ po.oheadr(constants.NOUT, "Numerics")
+ po.ocmmnt(constants.NOUT, "PROCESS has performed an evaluation run.")
+ po.oblnkl(constants.NOUT)
# Evaluate objective function
norm_objf = objective_function(numerics.minmax)
po.ovarre(
- constants.mfile, "Normalised objective function", "(norm_objf)", norm_objf
+ constants.MFILE, "Normalised objective function", "(norm_objf)", norm_objf
)
# Print the residuals of the constraint equations
@@ -78,12 +74,12 @@ def output_evaluation():
"Normalised residual": residual_error,
}
- po.write(constants.nout, tabulate(table_data, headers="keys"))
+ po.write(constants.NOUT, tabulate(table_data, headers="keys"))
for i in range(numerics.neqns):
constraint_id = numerics.icc[i]
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"{labels[i]} normalised residue",
f"(eq_con{constraint_id:03d})",
residual_error[i],
@@ -92,7 +88,7 @@ def output_evaluation():
for i in range(numerics.nineqns):
constraint_id = numerics.icc[numerics.neqns + i]
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"{labels[numerics.neqns + i]}",
f"(ineq_con{constraint_id:03d})",
residual_error[numerics.neqns + i],
diff --git a/process/fusion_reactions.py b/process/fusion_reactions.py
index efc03edc94..88720b0b9c 100644
--- a/process/fusion_reactions.py
+++ b/process/fusion_reactions.py
@@ -4,8 +4,8 @@
import numpy as np
from scipy import integrate
+from process import constants
from process.data_structure import physics_variables
-from process.fortran import constants
from process.plasma_profiles import PlasmaProfile
logger = logging.getLogger(__name__)
@@ -221,7 +221,7 @@ def dt_reaction(self) -> None:
self.sigmav_dt_average = sigmav
# Reaction energy in MegaJoules [MJ]
- reaction_energy = constants.d_t_energy / 1.0e6
+ reaction_energy = constants.D_T_ENERGY / 1.0e6
# Calculate the fusion power density produced [MW/m^3]
fusion_power_density = (
@@ -234,11 +234,11 @@ def dt_reaction(self) -> None:
# Power densities for different particles [MW/m^3]
# Alpha particle gets approximately 20% of the fusion power
alpha_power_density = (
- 1.0 - constants.dt_neutron_energy_fraction
+ 1.0 - constants.DT_NEUTRON_ENERGY_FRACTION
) * fusion_power_density
pden_non_alpha_charged_mw = 0.0
neutron_power_density = (
- constants.dt_neutron_energy_fraction * fusion_power_density
+ constants.DT_NEUTRON_ENERGY_FRACTION * fusion_power_density
)
# Calculate the fusion rate density [reactions/m^3/second]
@@ -290,7 +290,7 @@ def dhe3_reaction(self) -> None:
)
# Reaction energy in MegaJoules [MJ]
- reaction_energy = constants.d_helium_energy / 1.0e6
+ reaction_energy = constants.D_HELIUM_ENERGY / 1.0e6
# Calculate the fusion power density produced [MW/m^3]
fusion_power_density = (
@@ -303,10 +303,10 @@ def dhe3_reaction(self) -> None:
# Power densities for different particles [MW/m^3]
# Alpha particle gets approximately 20% of the fusion power
alpha_power_density = (
- 1.0 - constants.dhelium_proton_energy_fraction
+ 1.0 - constants.DHELIUM_PROTON_ENERGY_FRACTION
) * fusion_power_density
pden_non_alpha_charged_mw = (
- constants.dhelium_proton_energy_fraction * fusion_power_density
+ constants.DHELIUM_PROTON_ENERGY_FRACTION * fusion_power_density
)
neutron_power_density = 0.0
@@ -360,7 +360,7 @@ def dd_helion_reaction(self) -> None:
)
# Reaction energy in MegaJoules [MJ]
- reaction_energy = constants.dd_helium_energy / 1.0e6
+ reaction_energy = constants.DD_HELIUM_ENERGY / 1.0e6
# Calculate the fusion power density produced [MW/m^3]
# The power density is scaled by the branching ratio to simulate the different
@@ -377,10 +377,10 @@ def dd_helion_reaction(self) -> None:
# Neutron particle gets approximately 75% of the fusion power
alpha_power_density = 0.0
pden_non_alpha_charged_mw = (
- 1.0 - constants.dd_neutron_energy_fraction
+ 1.0 - constants.DD_NEUTRON_ENERGY_FRACTION
) * fusion_power_density
neutron_power_density = (
- constants.dd_neutron_energy_fraction * fusion_power_density
+ constants.DD_NEUTRON_ENERGY_FRACTION * fusion_power_density
)
# Calculate the fusion rate density [reactions/m^3/second]
@@ -433,7 +433,7 @@ def dd_triton_reaction(self) -> None:
)
# Reaction energy in MegaJoules [MJ]
- reaction_energy = constants.dd_triton_energy / 1.0e6
+ reaction_energy = constants.DD_TRITON_ENERGY / 1.0e6
# Calculate the fusion power density produced [MW/m^3]
# The power density is scaled by the branching ratio to simulate the different
@@ -741,8 +741,8 @@ def set_fusion_powers(
pden_neutron_total_mw = pden_plasma_neutron_mw + (
(
(
- constants.dt_neutron_energy_fraction
- / (1.0 - constants.dt_neutron_energy_fraction)
+ constants.DT_NEUTRON_ENERGY_FRACTION
+ / (1.0 - constants.DT_NEUTRON_ENERGY_FRACTION)
)
* p_beam_alpha_mw
)
@@ -864,8 +864,8 @@ def beam_fusion(
beam_slow_time = (
1.99e19
* (
- constants.m_deuteron_amu * (1.0 - f_beam_tritium)
- + (constants.m_triton_amu * f_beam_tritium)
+ constants.M_DEUTERON_AMU * (1.0 - f_beam_tritium)
+ + (constants.M_TRITON_AMU * f_beam_tritium)
)
* (ten**1.5 / dene)
/ ion_electron_coulomb_log
@@ -876,14 +876,14 @@ def beam_fusion(
# Taken from J.W Sheffield, “The physics of magnetic fusion reactors,”
critical_energy_deuterium = (
14.8
- * constants.m_deuteron_amu
+ * constants.M_DEUTERON_AMU
* ten
* zeffai ** (2 / 3)
* (ion_electron_coulomb_log + 4.0)
/ ion_electron_coulomb_log
)
critical_energy_tritium = critical_energy_deuterium * (
- constants.m_triton_amu / constants.m_deuteron_amu
+ constants.M_TRITON_AMU / constants.M_DEUTERON_AMU
)
# Deuterium and tritium ion densities
@@ -1006,7 +1006,7 @@ def beamcalc(
deuterium_beam_density = (
beam_current_deuterium
* characteristic_deuterium_beam_slow_time
- / (constants.electron_charge * vol_plasma)
+ / (constants.ELECTRON_CHARGE * vol_plasma)
)
# Ratio of beam energy to critical energy for tritium
@@ -1021,7 +1021,7 @@ def beamcalc(
tritium_beam_density = (
beam_current_tritium
* characteristic_tritium_beam_slow_time
- / (constants.electron_charge * vol_plasma)
+ / (constants.ELECTRON_CHARGE * vol_plasma)
)
hot_beam_density = deuterium_beam_density + tritium_beam_density
@@ -1030,43 +1030,43 @@ def beamcalc(
# Re-arrange kinetic energy equation to find speed. Non-relativistic.
deuterium_critical_energy_speed = np.sqrt(
2.0
- * constants.kiloelectron_volt
+ * constants.KILOELECTRON_VOLT
* critical_energy_deuterium
- / (constants.atomic_mass_unit * constants.m_deuteron_amu)
+ / (constants.ATOMIC_MASS_UNIT * constants.M_DEUTERON_AMU)
)
# Find the speed of the tritium particle when it has the critical energy.
# Re-arrange kinetic energy equation to find speed. Non-relativistic.
tritium_critical_energy_speed = np.sqrt(
2.0
- * constants.kiloelectron_volt
+ * constants.KILOELECTRON_VOLT
* critical_energy_tritium
- / (constants.atomic_mass_unit * constants.m_triton_amu)
+ / (constants.ATOMIC_MASS_UNIT * constants.M_TRITON_AMU)
)
# Source term representing the number of ions born per unit time per unit volume.
# D.Baiquan et.al. “Fast ion pressure in fusion plasma,” Nuclear Fusion and Plasma Physics,
# vol. 9, no. 3, pp. 136-141, 2022, Available: https://fti.neep.wisc.edu/fti.neep.wisc.edu/pdf/fdm718.pdf
- source_deuterium = beam_current_deuterium / (constants.electron_charge * vol_plasma)
+ source_deuterium = beam_current_deuterium / (constants.ELECTRON_CHARGE * vol_plasma)
- source_tritium = beam_current_tritium / (constants.electron_charge * vol_plasma)
+ source_tritium = beam_current_tritium / (constants.ELECTRON_CHARGE * vol_plasma)
pressure_coeff_deuterium = (
- constants.m_deuteron_amu
- * constants.atomic_mass_unit
+ constants.M_DEUTERON_AMU
+ * constants.ATOMIC_MASS_UNIT
* beam_slow_time
* deuterium_critical_energy_speed**2
* source_deuterium
- / (constants.kiloelectron_volt * 3.0)
+ / (constants.KILOELECTRON_VOLT * 3.0)
)
pressure_coeff_tritium = (
- constants.m_triton_amu
- * constants.atomic_mass_unit
+ constants.M_TRITON_AMU
+ * constants.ATOMIC_MASS_UNIT
* beam_slow_time
* tritium_critical_energy_speed**2
* source_tritium
- / (constants.kiloelectron_volt * 3.0)
+ / (constants.KILOELECTRON_VOLT * 3.0)
)
# Fast Ion Pressure
@@ -1089,11 +1089,11 @@ def beamcalc(
) / hot_beam_density
hot_deuterium_rate = 1e-4 * beam_reaction_rate(
- constants.m_deuteron_amu, deuterium_critical_energy_speed, e_beam_kev
+ constants.M_DEUTERON_AMU, deuterium_critical_energy_speed, e_beam_kev
)
hot_tritium_rate = 1e-4 * beam_reaction_rate(
- constants.m_triton_amu, tritium_critical_energy_speed, e_beam_kev
+ constants.M_TRITON_AMU, tritium_critical_energy_speed, e_beam_kev
)
deuterium_beam_alpha_power = alpha_power_beam(
@@ -1200,7 +1200,7 @@ def alpha_power_beam(
beam_ion_desnity
* plasma_ion_desnity
* sigv
- * (constants.dt_alpha_energy / 1e6)
+ * (constants.DT_ALPHA_ENERGY / 1e6)
* vol_plasma
* ratio
)
@@ -1236,9 +1236,9 @@ def beam_reaction_rate(
# Find the speed of the beam particle when it has the critical energy.
# Re-arrange kinetic energy equation to find speed. Non-relativistic.
beam_velocity = np.sqrt(
- (beam_energy_keV * constants.kiloelectron_volt)
+ (beam_energy_keV * constants.KILOELECTRON_VOLT)
* 2.0
- / (relative_mass_ion * constants.atomic_mass_unit)
+ / (relative_mass_ion * constants.ATOMIC_MASS_UNIT)
)
relative_velocity = beam_velocity / critical_velocity
@@ -1288,7 +1288,7 @@ def _hot_beam_fusion_reaction_rate_integrand(
beam_velcoity = critical_velocity * velocity_ratio
# Calculate the beam kinetic energy per amu and normalise to keV
- xvcs = beam_velcoity**2 * constants.atomic_mass_unit / (constants.kiloelectron_volt)
+ xvcs = beam_velcoity**2 * constants.ATOMIC_MASS_UNIT / (constants.KILOELECTRON_VOLT)
# Calculate the fusion reaction cross-section from beam kinetic energy
cross_section = _beam_fusion_cross_section(xvcs)
@@ -1328,7 +1328,7 @@ def _beam_fusion_cross_section(vrelsq: float) -> float:
a5 = 4.09e2
# Beam kinetic energy
- e_beam_kev = 0.5 * constants.m_deuteron_amu * vrelsq
+ e_beam_kev = 0.5 * constants.M_DEUTERON_AMU * vrelsq
# Set limits on cross-section at low and high beam energies
if e_beam_kev < 10.0:
diff --git a/process/fw.py b/process/fw.py
index 1f3453a684..2be2d6e36a 100644
--- a/process/fw.py
+++ b/process/fw.py
@@ -2,17 +2,17 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.coolprop_interface import FluidProperties
from process.data_structure import blanket_library, build_variables, fwbs_variables
-from process.fortran import constants
logger = logging.getLogger(__name__)
class Fw:
def __init__(self) -> None:
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def run(self):
(
diff --git a/process/hcpb.py b/process/hcpb.py
index fbeae9e4ad..cf7a648e4e 100644
--- a/process/hcpb.py
+++ b/process/hcpb.py
@@ -3,6 +3,7 @@
import numpy as np
import process.blanket_library as blanket_library
+from process import constants
from process import (
process_output as po,
)
@@ -21,7 +22,6 @@
tfcoil_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import constants
logger = logging.getLogger(__name__)
@@ -349,7 +349,7 @@ def component_masses(self):
# First wall armour mass (kg)
fwbs_variables.fw_armour_mass = (
- fwbs_variables.fw_armour_vol * constants.den_tungsten
+ fwbs_variables.fw_armour_vol * constants.DEN_TUNGSTEN
)
if fwbs_variables.breeder_f < 1.0e-10:
@@ -444,7 +444,7 @@ def nuclear_heating_magnets(self, output: bool):
# Calculate smeared densities of blanket sections
# gaseous He coolant in armour, FW & blanket: He mass is neglected
- ccfe_hcpb_module.armour_density = constants.den_tungsten * (1.0 - vffwm)
+ ccfe_hcpb_module.armour_density = constants.DEN_TUNGSTEN * (1.0 - vffwm)
ccfe_hcpb_module.fw_density = fwbs_variables.den_steel * (1.0 - vffwm)
ccfe_hcpb_module.blanket_density = (
fwbs_variables.m_blkt_total / fwbs_variables.vol_blkt_total
diff --git a/process/ife.py b/process/ife.py
index 2287fda892..169d585701 100644
--- a/process/ife.py
+++ b/process/ife.py
@@ -7,7 +7,7 @@
import numpy as np
-from process import process_output
+from process import constants, process_output
from process.data_structure import (
build_variables,
buildings_variables,
@@ -20,7 +20,6 @@
vacuum_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import constants
MATERIALS = [
"void",
@@ -52,7 +51,7 @@ def __init__(self, availability, costs) -> None:
:type costs: process.costs.Costs
"""
- self.outfile: int = constants.nout
+ self.outfile: int = constants.NOUT
self.availability = availability
self.costs = costs
@@ -892,7 +891,7 @@ def bld2019(self):
# Velocity
vel = np.sqrt(
2.0
- * constants.acceleration_gravity
+ * constants.ACCELERATION_GRAVITY
* (ife_variables.chdzu + ife_variables.bldzu)
)
@@ -917,7 +916,7 @@ def bld2019(self):
ife_variables.lipmw = (
1e-6
* mdot
- * constants.acceleration_gravity
+ * constants.ACCELERATION_GRAVITY
* (
ife_variables.chdzl
+ ife_variables.chdzu
@@ -1980,7 +1979,7 @@ def ifepw1(self):
1.0e6
* heat_transport_variables.p_cryo_plant_electric_mw
* (0.13 * 4.5)
- / (constants.temp_room - 4.5)
+ / (constants.TEMP_ROOM - 4.5)
)
def ifepw2(self, output: bool = False):
diff --git a/process/impurity_radiation.py b/process/impurity_radiation.py
index 20d1f1f344..c6084ca06b 100644
--- a/process/impurity_radiation.py
+++ b/process/impurity_radiation.py
@@ -7,9 +7,9 @@
import numpy as np
from scipy import integrate
+from process import constants
from process.data_structure import impurity_radiation_module
from process.exceptions import ProcessError, ProcessValueError
-from process.fortran import constants
logger = logging.getLogger(__name__)
@@ -34,7 +34,7 @@ def initialise_imprad():
no=1,
label=impurity_radiation_module.imp_label[0],
z=1,
- amass=constants.m_protium_amu, # 1.00782503223 1H
+ amass=constants.M_PROTIUM_AMU, # 1.00782503223 1H
frac=frac,
len_tab=table_length,
error=errorflag,
@@ -47,7 +47,7 @@ def initialise_imprad():
2,
impurity_radiation_module.imp_label[1],
2,
- constants.m_helium_amu, # 4.002602 (3He,4He) Average mass
+ constants.M_HELIUM_AMU, # 4.002602 (3He,4He) Average mass
frac,
table_length,
errorflag,
@@ -58,7 +58,7 @@ def initialise_imprad():
3,
impurity_radiation_module.imp_label[2],
4,
- constants.m_beryllium_amu, # 9.0121831 9Be
+ constants.M_BERYLLIUM_AMU, # 9.0121831 9Be
frac,
table_length,
errorflag,
@@ -69,7 +69,7 @@ def initialise_imprad():
4,
impurity_radiation_module.imp_label[3],
6,
- constants.m_carbon_amu, # 12.0096, (12C,13C,14C) Average mass
+ constants.M_CARBON_AMU, # 12.0096, (12C,13C,14C) Average mass
frac,
table_length,
errorflag,
@@ -80,7 +80,7 @@ def initialise_imprad():
5,
impurity_radiation_module.imp_label[4],
7,
- constants.m_nitrogen_amu, # 14.00643, (14N,15N) Average mass
+ constants.M_NITROGEN_AMU, # 14.00643, (14N,15N) Average mass
frac,
table_length,
errorflag,
@@ -91,7 +91,7 @@ def initialise_imprad():
6,
impurity_radiation_module.imp_label[5],
8,
- constants.m_oxygen_amu, # 15.99903, (16O,17O,18O) Average mass
+ constants.M_OXYGEN_AMU, # 15.99903, (16O,17O,18O) Average mass
frac,
table_length,
errorflag,
@@ -102,7 +102,7 @@ def initialise_imprad():
7,
impurity_radiation_module.imp_label[6],
10,
- constants.m_neon_amu, # 20.1797 (20Ne,21Ne,22Ne) Average mass
+ constants.M_NEON_AMU, # 20.1797 (20Ne,21Ne,22Ne) Average mass
frac,
table_length,
errorflag,
@@ -113,7 +113,7 @@ def initialise_imprad():
8,
impurity_radiation_module.imp_label[7],
14,
- constants.m_silicon_amu, # 28.084 (28Si,29Si,30Si) Average mass
+ constants.M_SILICON_AMU, # 28.084 (28Si,29Si,30Si) Average mass
frac,
table_length,
errorflag,
@@ -124,7 +124,7 @@ def initialise_imprad():
9,
impurity_radiation_module.imp_label[8],
18,
- constants.m_argon_amu, # 39.948 (40Ar,36Ar,38Ar) Average mass
+ constants.M_ARGON_AMU, # 39.948 (40Ar,36Ar,38Ar) Average mass
frac,
table_length,
errorflag,
@@ -135,7 +135,7 @@ def initialise_imprad():
10,
impurity_radiation_module.imp_label[9],
26,
- constants.m_iron_amu, # 55.845 (56Fe,54Fe,57Fe,58Fe) Average mass
+ constants.M_IRON_AMU, # 55.845 (56Fe,54Fe,57Fe,58Fe) Average mass
frac,
table_length,
errorflag,
@@ -146,7 +146,7 @@ def initialise_imprad():
11,
impurity_radiation_module.imp_label[10],
28,
- constants.m_nickel_amu, # 58.6934 (58Ni,60Ni,61Ni,62Ni,64Ni) Average mass
+ constants.M_NICKEL_AMU, # 58.6934 (58Ni,60Ni,61Ni,62Ni,64Ni) Average mass
frac,
table_length,
errorflag,
@@ -157,7 +157,7 @@ def initialise_imprad():
12,
impurity_radiation_module.imp_label[11],
36,
- constants.m_krypton_amu, # 83.798 (84Kr,86Kr,82Kr,80Kr,78Kr) Average mass
+ constants.M_KRYPTON_AMU, # 83.798 (84Kr,86Kr,82Kr,80Kr,78Kr) Average mass
frac,
table_length,
errorflag,
@@ -168,7 +168,7 @@ def initialise_imprad():
13,
impurity_radiation_module.imp_label[12],
54,
- constants.m_xenon_amu, # 131.293 (132Xe,129Xe,131Xe,134Xe,136Xe) Average mass
+ constants.M_XENON_AMU, # 131.293 (132Xe,129Xe,131Xe,134Xe,136Xe) Average mass
frac,
table_length,
errorflag,
@@ -179,7 +179,7 @@ def initialise_imprad():
14,
impurity_radiation_module.imp_label[13],
74,
- constants.m_tungsten_amu, # 183.84 (184W,186W,182W,183W,180W) Average mass
+ constants.M_TUNGSTEN_AMU, # 183.84 (184W,186W,182W,183W,180W) Average mass
frac,
table_length,
errorflag,
diff --git a/process/init.py b/process/init.py
index b506f01520..e9f7d22679 100644
--- a/process/init.py
+++ b/process/init.py
@@ -9,7 +9,7 @@
import process.fortran as fortran
import process.iteration_variables as iteration_variables
import process.process_output as process_output
-from process import data_structure
+from process import constants, data_structure
from process.constraints import ConstraintManager
from process.data_structure.blanket_library import init_blanket_library
from process.data_structure.build_variables import init_build_variables
@@ -75,9 +75,8 @@ def init_process():
# Initialise the program variables
iteration_variables.initialise_iteration_variables()
- # Initialise the Fortran file specifiers
- # (creating and opening the files in the process)
- fortran.init_module.open_files()
+ # Creating and open the files MFile and OUTFile
+ process_output.OutputFileManager.open_files()
# Input any desired new initial values
inputs = parse_input_file()
@@ -133,7 +132,7 @@ def get_git_summary() -> tuple[str, str]:
def run_summary():
"""Write a summary of the PROCESS run to the output file and MFile"""
# Outfile and terminal #
- for outfile in [fortran.constants.nout, fortran.constants.iotty]:
+ for outfile in [constants.NOUT, constants.IOTTY]:
# PROCESS code header
process_output.oblnkl(outfile)
process_output.ostars(outfile, 110)
@@ -166,12 +165,12 @@ def run_summary():
process_output.ocmmnt(outfile, f"Computer : {machine}")
process_output.ocmmnt(outfile, f"Directory : {Path.cwd()}")
- fileprefix = f2py_compatible_to_string(fortran.global_variables.fileprefix)
+ fileprefix = data_structure.global_variables.fileprefix
process_output.ocmmnt(
outfile,
f"Input : {fileprefix}",
)
- runtitle = f2py_compatible_to_string(fortran.global_variables.runtitle)
+ runtitle = data_structure.global_variables.runtitle
process_output.ocmmnt(
outfile,
f"Run title : {runtitle}",
@@ -179,7 +178,7 @@ def run_summary():
process_output.ocmmnt(
outfile,
- f"Run type : Reactor concept design: {f2py_compatible_to_string(fortran.global_variables.icase)}, (c) UK Atomic Energy Authority",
+ f"Run type : Reactor concept design: {data_structure.global_variables.icase}, (c) UK Atomic Energy Authority",
)
process_output.oblnkl(outfile)
@@ -202,7 +201,8 @@ def run_summary():
# If optimising, write objective function and convergence parameter
if fortran.numerics.ioptimz == 1:
process_output.ocmmnt(
- outfile, f"Max iterations : {fortran.global_variables.maxcal.item()}"
+ outfile,
+ f"Max iterations : {data_structure.global_variables.maxcal}",
)
if fortran.numerics.minmax > 0:
@@ -228,7 +228,7 @@ def run_summary():
process_output.ostars(outfile, 110)
# MFile #
- mfile = fortran.constants.mfile
+ mfile = constants.MFILE
process_output.ovarst(mfile, "PROCESS version", "(procver)", f'"{version}"')
process_output.ovarst(mfile, "Date of run", "(date)", f'"{date_string}"')
@@ -264,7 +264,7 @@ def init_all_module_vars():
init_cost_variables()
init_divertor_variables()
init_fwbs_variables()
- fortran.global_variables.init_global_variables()
+ data_structure.global_variables.init_global_variables()
init_ccfe_hcpb_module()
init_heat_transport_variables()
init_ife_variables()
@@ -277,7 +277,7 @@ def init_all_module_vars():
init_stellarator_variables()
init_tfcoil_variables()
init_times_variables()
- fortran.constants.init_constants()
+ constants.init_constants()
init_current_drive_variables()
init_primary_pumping_variables()
init_pfcoil_variables()
@@ -585,7 +585,7 @@ def check_process(inputs): # noqa: ARG001
# Tight aspect ratio options (ST)
if data_structure.physics_variables.itart == 1:
- fortran.global_variables.icase = "Tight aspect ratio tokamak model"
+ data_structure.global_variables.icase = "Tight aspect ratio tokamak model"
# Disabled Forcing that no inboard breeding blanket is used
# Disabled i_blkt_inboard = 0
@@ -778,7 +778,7 @@ def check_process(inputs): # noqa: ARG001
# Pulsed power plant model
if data_structure.pulse_variables.i_pulsed_plant == 1:
- fortran.global_variables.icase = "Pulsed tokamak model"
+ data_structure.global_variables.icase = "Pulsed tokamak model"
else:
data_structure.buildings_variables.esbldgm3 = 0.0
@@ -1206,6 +1206,6 @@ def set_active_constraints():
def set_device_type():
if data_structure.ife_variables.ife == 1:
- fortran.global_variables.icase = "Inertial Fusion model"
+ data_structure.global_variables.icase = "Inertial Fusion model"
elif data_structure.stellarator_variables.istell != 0:
- fortran.global_variables.icase = "Stellarator model"
+ data_structure.global_variables.icase = "Stellarator model"
diff --git a/process/input.py b/process/input.py
index 8e717dc01a..fe3a5567c6 100644
--- a/process/input.py
+++ b/process/input.py
@@ -8,14 +8,11 @@
from typing import Any
from warnings import warn
+import process
import process.data_structure as data_structure
import process.fortran as fortran
from process.constraints import ConstraintManager
from process.exceptions import ProcessValidationError, ProcessValueError
-from process.utilities.f2py_string_patch import (
- f2py_compatible_to_string,
- string_to_f2py_compatible,
-)
NumberType = int | float
ValidInputTypes = NumberType | str
@@ -95,15 +92,15 @@ def __post_init__(self):
INPUT_VARIABLES = {
- "runtitle": InputVariable(fortran.global_variables, str),
- "verbose": InputVariable(fortran.global_variables, int, choices=[0, 1]),
- "run_tests": InputVariable(fortran.global_variables, int, choices=[0, 1]),
+ "runtitle": InputVariable(data_structure.global_variables, str),
+ "verbose": InputVariable(data_structure.global_variables, int, choices=[0, 1]),
+ "run_tests": InputVariable(data_structure.global_variables, int, choices=[0, 1]),
"ioptimz": InputVariable(fortran.numerics, int, choices=[1, -2]),
"epsvmc": InputVariable(fortran.numerics, float, range=(0.0, 1.0)),
"boundl": InputVariable(fortran.numerics, float, array=True),
"boundu": InputVariable(fortran.numerics, float, array=True),
"epsfcn": InputVariable(fortran.numerics, float, range=(0.0, 1.0)),
- "maxcal": InputVariable(fortran.global_variables, int, range=(0, 10000)),
+ "maxcal": InputVariable(data_structure.global_variables, int, range=(0, 10000)),
"minmax": InputVariable(fortran.numerics, int),
"neqns": InputVariable(
fortran.numerics, int, range=(1, ConstraintManager.num_constraints())
@@ -452,7 +449,7 @@ def __post_init__(self):
"dz_vv_upper": InputVariable(
data_structure.build_variables, float, range=(0.0, 10.0)
),
- "den_aluminium": InputVariable(fortran.constants, float, range=(2500.0, 30000.0)),
+ "den_aluminium": InputVariable(process.constants, float, range=(2500.0, 30000.0)),
"den_tf_coil_case": InputVariable(
data_structure.tfcoil_variables, float, range=(1000.0, 100000.0)
),
@@ -462,7 +459,7 @@ def __post_init__(self):
"den_tf_wp_turn_insulation": InputVariable(
data_structure.tfcoil_variables, float, range=(500.0, 10000.0)
),
- "den_copper": InputVariable(fortran.constants, float, range=(8000.0, 10000.0)),
+ "den_copper": InputVariable(process.constants, float, range=(8000.0, 10000.0)),
"declblkt": InputVariable(data_structure.fwbs_variables, float, range=(0.01, 0.2)),
"declfw": InputVariable(data_structure.fwbs_variables, float, range=(0.01, 0.2)),
"declshld": InputVariable(data_structure.fwbs_variables, float, range=(0.01, 0.2)),
@@ -2310,7 +2307,7 @@ def __post_init__(self):
def parse_input_file():
- input_file = f2py_compatible_to_string(fortran.global_variables.fileprefix)
+ input_file = data_structure.global_variables.fileprefix
input_file_path = Path("IN.DAT")
if input_file != "":
@@ -2501,11 +2498,7 @@ def set_scalar_variable(name: str, value: ValidInputTypes, config: InputVariable
)
raise ProcessValueError(error_msg)
- # Can remove once global_variables is in Python
- if name == "runtitle":
- setattr(config.module, name, string_to_f2py_compatible(current_value, value))
- else:
- setattr(config.module, name, value)
+ setattr(config.module, name, value)
def set_array_variable(name: str, value: str, array_index: int, config: InputVariable):
diff --git a/process/iteration_variables.py b/process/iteration_variables.py
index 8616306bc4..3fba704434 100644
--- a/process/iteration_variables.py
+++ b/process/iteration_variables.py
@@ -9,7 +9,6 @@
import process.fortran as fortran
from process.exceptions import ProcessValueError
from process.utilities.f2py_string_patch import (
- f2py_compatible_to_string,
string_to_f2py_compatible,
)
@@ -469,8 +468,8 @@ def load_iteration_variables():
# warn of the iteration variable is also a scan variable because this will cause
# the optimiser and scan to overwrite the same variable and conflict
if iteration_variable.name in (
- f2py_compatible_to_string(fortran.global_variables.vlabel),
- f2py_compatible_to_string(fortran.global_variables.vlabel_2),
+ data_structure.global_variables.vlabel,
+ data_structure.global_variables.vlabel_2,
):
warn(
(
diff --git a/process/log.py b/process/log.py
index 6dad5532ae..f4368806db 100644
--- a/process/log.py
+++ b/process/log.py
@@ -9,7 +9,7 @@
from logging import Handler
import process.process_output as process_output
-from process.fortran import constants
+from process import constants
class ProcessLogHandler(Handler):
@@ -74,7 +74,7 @@ def show_errors(file_unit: int):
print(warning_string)
process_output.write(file_unit, warning_string)
process_output.ovarre(
- constants.mfile,
+ constants.MFILE,
"Error status",
"(error_status)",
0 if logging_model_handler.num_logs() == 0 else 2,
diff --git a/process/main.py b/process/main.py
index e769006b43..1435abdad5 100644
--- a/process/main.py
+++ b/process/main.py
@@ -51,6 +51,7 @@
import process.data_structure as data_structure
import process.fortran as fortran
import process.init as init
+from process import constants
from process.availability import Availability
from process.blanket_library import BlanketLibrary
from process.build import Build
@@ -93,6 +94,7 @@
from process.plasma_geometry import PlasmaGeom
from process.plasma_profiles import PlasmaProfile
from process.power import Power
+from process.process_output import OutputFileManager, oheadr
from process.pulse import Pulse
from process.resistive_tf_coil import AluminiumTFCoil, CopperTFCoil, ResistiveTFCoil
from process.scan import Scan
@@ -100,7 +102,6 @@
from process.structure import Structure
from process.superconducting_tf_coil import SuperconductingTFCoil
from process.tf_coil import TFCoil
-from process.utilities.f2py_string_patch import string_to_f2py_compatible
from process.vacuum import Vacuum
from process.water_use import WaterUse
@@ -422,11 +423,7 @@ def set_input(self):
)
# Set the input file in the Fortran
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- str(self.input_path.resolve()),
- except_length=True,
- )
+ data_structure.global_variables.fileprefix = str(self.input_path.resolve())
def set_output(self):
"""Set the output file name.
@@ -434,9 +431,7 @@ def set_output(self):
Set Path object on the Process object, and set the prefix in the Fortran.
"""
self.output_path = Path(self.filename_prefix + "OUT.DAT")
- fortran.global_variables.output_prefix = string_to_f2py_compatible(
- fortran.global_variables.output_prefix, self.filename_prefix
- )
+ data_structure.global_variables.output_prefix = self.filename_prefix
def set_mfile(self):
"""Set the mfile filename."""
@@ -474,7 +469,7 @@ def run_scan(self):
def show_errors(self):
"""Report all informational/error messages encountered."""
- show_errors(fortran.constants.nout)
+ show_errors(constants.NOUT)
def finish(self):
"""Run the finish subroutine to close files open in the Fortran.
@@ -482,7 +477,10 @@ def finish(self):
Files being handled by Fortran must be closed before attempting to
write to them using Python, otherwise only parts are written.
"""
- fortran.init_module.finish()
+ oheadr(constants.NOUT, "End of PROCESS Output")
+ oheadr(constants.IOTTY, "End of PROCESS Output")
+ oheadr(constants.NOUT, "Copy of PROCESS Input Follows")
+ OutputFileManager.finish()
def append_input(self):
"""Append the input file to the output file and mfile."""
diff --git a/process/pfcoil.py b/process/pfcoil.py
index fb1c2eb106..e471d650db 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -8,7 +8,7 @@
from scipy.special import ellipe, ellipk
import process.superconductors as superconductors
-from process import fortran as ft
+from process import constants
from process import process_output as op
from process.data_structure import build_variables as bv
from process.data_structure import constraint_variables as ctv
@@ -20,20 +20,18 @@
from process.data_structure import tfcoil_variables as tfv
from process.data_structure import times_variables as tv
from process.exceptions import ProcessValueError
-from process.fortran import constants, numerics
+from process.fortran import numerics
logger = logging.getLogger(__name__)
-RMU0 = ft.constants.rmu0
-
class PFCoil:
"""Calculate poloidal field coil system parameters."""
def __init__(self, cs_fatigue) -> None:
"""Initialise Fortran module variables."""
- self.outfile = ft.constants.nout # output file unit
- self.mfile = ft.constants.mfile # mfile file unit
+ self.outfile = constants.NOUT # output file unit
+ self.mfile = constants.MFILE # mfile file unit
pfcoil_variables.init_pfcoil_module()
self.cs_fatigue = cs_fatigue
@@ -556,8 +554,8 @@ def pfcoil(self):
ddics = (
4.0e-7
- * constants.pi
- * constants.pi
+ * constants.PI
+ * constants.PI
* (
(bv.dr_bore * bv.dr_bore)
+ (bv.dr_cs * bv.dr_cs) / 6.0e0
@@ -806,7 +804,7 @@ def pfcoil(self):
rll = (
2.0e0
- * constants.pi
+ * constants.PI
* pfcoil_variables.r_pf_coil_middle[i]
* pfcoil_variables.n_pf_coil_turns[i]
)
@@ -904,7 +902,7 @@ def pfcoil(self):
pfcoil_variables.m_pf_coil_structure[i] = (
areaspf
* 2.0e0
- * constants.pi
+ * constants.PI
* pfcoil_variables.r_pf_coil_middle[i]
* fwbsv.den_steel
)
@@ -1326,22 +1324,22 @@ def ohcalc(self):
# CS coil turn geometry calculation - stadium shape
# Literature: https://doi.org/10.1016/j.fusengdes.2017.04.052
pfcoil_variables.radius_cs_turn_cable_space = -(
- (pfcoil_variables.dr_cs_turn - pfcoil_variables.dz_cs_turn) / constants.pi
+ (pfcoil_variables.dr_cs_turn - pfcoil_variables.dz_cs_turn) / constants.PI
) + math.sqrt(
(
(
(pfcoil_variables.dr_cs_turn - pfcoil_variables.dz_cs_turn)
- / constants.pi
+ / constants.PI
)
** 2
)
+ (
(
(pfcoil_variables.dr_cs_turn * pfcoil_variables.dz_cs_turn)
- - (4 - constants.pi) * (pfcoil_variables.r_out_cst**2)
+ - (4 - constants.PI) * (pfcoil_variables.r_out_cst**2)
- (pfcoil_variables.a_cs_turn * pfcoil_variables.f_a_cs_steel)
)
- / constants.pi
+ / constants.PI
)
)
@@ -1466,7 +1464,7 @@ def ohcalc(self):
pfcoil_variables.m_pf_coil_structure[pfcoil_variables.n_cs_pf_coils - 1] = (
areaspf
* 2.0e0
- * constants.pi
+ * constants.PI
* pfcoil_variables.r_pf_coil_middle[pfcoil_variables.n_cs_pf_coils - 1]
* fwbsv.den_steel
)
@@ -1487,7 +1485,7 @@ def ohcalc(self):
pfcoil_variables.awpoh
* (1.0e0 - pfcoil_variables.f_a_cs_void)
* 2.0e0
- * constants.pi
+ * constants.PI
* pfcoil_variables.r_pf_coil_middle[pfcoil_variables.n_cs_pf_coils - 1]
* tfv.dcond[pfcoil_variables.i_cs_superconductor - 1]
)
@@ -1496,7 +1494,7 @@ def ohcalc(self):
pfcoil_variables.awpoh
* (1.0e0 - pfcoil_variables.f_a_cs_void)
* 2.0e0
- * constants.pi
+ * constants.PI
* pfcoil_variables.r_pf_coil_middle[pfcoil_variables.n_cs_pf_coils - 1]
* constants.den_copper
)
@@ -1588,7 +1586,7 @@ def ohcalc(self):
pfcoil_variables.p_cs_resistive_flat_top = (
2.0e0
- * constants.pi
+ * constants.PI
* pfcoil_variables.r_cs_middle
* pfcoil_variables.rho_pf_coil
/ (
@@ -1827,7 +1825,7 @@ def bfmax(self, rj, a, b, h):
# Fits are for 1 < alpha < 2 , and 0.5 < beta < very large
b0 = (
rj
- * constants.rmu0
+ * constants.RMU0
* h
* math.log(
(alpha + math.sqrt(alpha**2 + beta**2))
@@ -1836,7 +1834,7 @@ def bfmax(self, rj, a, b, h):
)
if beta > 3.0:
- b1 = constants.rmu0 * rj * (b - a)
+ b1 = constants.RMU0 * rj * (b - a)
f = (3.0 / beta) ** 2
bfmax = f * b0 * (1.007 + (alpha - 1.0) * 0.0055) + (1.0 - f) * b1
@@ -2070,7 +2068,7 @@ def axial_stress(self):
k2b2 = (4.0e0 * b**2) / (4.0e0 * b**2 + 4.0e0 * hl**2)
# term 1
- axial_term_1 = -(constants.rmu0 / 2.0e0) * (ni / (2.0e0 * hl)) ** 2
+ axial_term_1 = -(constants.RMU0 / 2.0e0) * (ni / (2.0e0 * hl)) ** 2
# term 2
ekb2_1 = ellipk(kb2)
@@ -2090,7 +2088,7 @@ def axial_stress(self):
axial_force = axial_term_1 * (axial_term_2 - axial_term_3)
# axial area [m2]
- area_ax = constants.pi * (
+ area_ax = constants.PI * (
pfcoil_variables.r_pf_coil_outer[pfcoil_variables.n_cs_pf_coils - 1] ** 2
- pfcoil_variables.r_pf_coil_inner[pfcoil_variables.n_cs_pf_coils - 1] ** 2
)
@@ -2339,9 +2337,9 @@ def induct(self, output):
rl = abs(
pfcoil_variables.z_pf_coil_upper[k]
- pfcoil_variables.z_pf_coil_lower[k]
- ) / math.sqrt(constants.pi)
+ ) / math.sqrt(constants.PI)
pfcoil_variables.ind_pf_cs_plasma_mutual[k, k] = (
- constants.rmu0
+ constants.RMU0
* pfcoil_variables.n_pf_coil_turns[k] ** 2
* pfcoil_variables.r_pf_coil_middle[k]
* (
@@ -3708,19 +3706,19 @@ def bfield(rc, zc, cc, rp, zp):
# Mutual inductances
- xc[i] = 0.5 * RMU0 * sd * ((2.0 - s) * xk - 2.0 * xe)
+ xc[i] = 0.5 * constants.RMU0 * sd * ((2.0 - s) * xk - 2.0 * xe)
# Radial, vertical fields
brx = (
- RMU0
+ constants.RMU0
* cc[i]
* dz
/ (2 * np.pi * rp * sd)
* (-xk + (rc[i] ** 2 + rp**2 + zs) / (dr**2 + zs) * xe)
)
bzx = (
- RMU0
+ constants.RMU0
* cc[i]
/ (2 * np.pi * sd)
* (xk + (rc[i] ** 2 - rp**2 - zs) / (dr**2 + zs) * xe)
diff --git a/process/physics.py b/process/physics.py
index 6b527fc72d..c07c822f50 100644
--- a/process/physics.py
+++ b/process/physics.py
@@ -12,9 +12,8 @@
import process.impurity_radiation as impurity_radiation
import process.l_h_transition as transition
import process.physics_functions as physics_funcs
-from process import (
- process_output as po,
-)
+from process import constants
+from process import process_output as po
from process.data_structure import (
build_variables,
constraint_variables,
@@ -29,13 +28,10 @@
times_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import constants, numerics
+from process.fortran import numerics
logger = logging.getLogger(__name__)
-ELECTRON_CHARGE: float = constants.electron_charge.item()
-RMU0: float = constants.rmu0.item()
-
@nb.jit(nopython=True, cache=True)
def rether(alphan, alphat, dene, dlamie, te, ti, zeffai):
@@ -130,7 +126,7 @@ def calculate_volt_second_requirements(
"""
# Plasma internal inductance
- ind_plasma_internal = RMU0 * rmajor * ind_plasma_internal_norm / 2.0
+ ind_plasma_internal = constants.RMU0 * rmajor * ind_plasma_internal_norm / 2.0
# Internal plasma flux (V-s) component
vs_plasma_internal = ind_plasma_internal * plasma_current
@@ -138,7 +134,7 @@ def calculate_volt_second_requirements(
# Start-up resistive component
# Uses ITER formula without the 10 V-s add-on
- vs_res_ramp = ejima_coeff * RMU0 * plasma_current * rmajor
+ vs_res_ramp = ejima_coeff * constants.RMU0 * plasma_current * rmajor
# ======================================================================
@@ -150,7 +146,7 @@ def calculate_volt_second_requirements(
beps = 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6)
ind_plasma_external = (
- rmajor * RMU0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa)
+ rmajor * constants.RMU0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa)
)
# ======================================================================
@@ -688,7 +684,9 @@ def _nevins_integral(
"""
# Compute average electron beta
- betae = dene * te * 1.0e3 * ELECTRON_CHARGE / (bt**2 / (2.0 * RMU0))
+ betae = (
+ dene * te * 1.0e3 * constants.ELECTRON_CHARGE / (bt**2 / (2.0 * constants.RMU0))
+ )
nabla = rminor * np.sqrt(y) / rmajor
x = (1.46 * np.sqrt(nabla) + 2.4 * nabla) / (1.0 - nabla) ** 1.5
@@ -1529,7 +1527,7 @@ def _trapped_particle_fraction_sauter(
class Physics:
def __init__(self, plasma_profile, current_drive):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
self.plasma_profile = plasma_profile
self.current_drive = current_drive
@@ -1756,7 +1754,7 @@ def physics(self):
* physics_variables.beta_thermal
* physics_variables.btot
* physics_variables.btot
- / (2.0e0 * constants.rmu0)
+ / (2.0e0 * constants.RMU0)
* physics_variables.vol_plasma
)
@@ -1766,7 +1764,7 @@ def physics(self):
* physics_variables.beta
* physics_variables.btot
* physics_variables.btot
- / (2.0e0 * constants.rmu0)
+ / (2.0e0 * constants.RMU0)
* physics_variables.vol_plasma
)
@@ -2159,27 +2157,27 @@ def physics(self):
physics_variables.fusden_plasma
+ 1.0e6
* physics_variables.p_beam_alpha_mw
- / (constants.dt_alpha_energy)
+ / (constants.DT_ALPHA_ENERGY)
/ physics_variables.vol_plasma
)
physics_variables.fusden_alpha_total = (
physics_variables.fusden_plasma_alpha
+ 1.0e6
* physics_variables.p_beam_alpha_mw
- / (constants.dt_alpha_energy)
+ / (constants.DT_ALPHA_ENERGY)
/ physics_variables.vol_plasma
)
physics_variables.p_dt_total_mw = (
physics_variables.p_plasma_dt_mw
- + (1.0 / (1.0 - constants.dt_neutron_energy_fraction))
+ + (1.0 / (1.0 - constants.DT_NEUTRON_ENERGY_FRACTION))
* physics_variables.p_beam_alpha_mw
)
physics_variables.p_beam_neutron_mw = physics_variables.p_beam_alpha_mw * (
- constants.dt_neutron_energy_fraction
- / (1 - constants.dt_neutron_energy_fraction)
+ constants.DT_NEUTRON_ENERGY_FRACTION
+ / (1 - constants.DT_NEUTRON_ENERGY_FRACTION)
)
physics_variables.p_beam_dt_mw = physics_variables.p_beam_alpha_mw * (
- 1 / (1 - constants.dt_neutron_energy_fraction)
+ 1 / (1 - constants.DT_NEUTRON_ENERGY_FRACTION)
)
else:
# If no beams present then the total alpha rates and power are the same as the plasma values
@@ -3305,9 +3303,9 @@ def plasma_composition() -> None:
# Average atomic masses of injected fuel species
physics_variables.m_fuel_amu = (
- (constants.m_deuteron_amu * physics_variables.f_deuterium)
- + (constants.m_triton_amu * physics_variables.f_tritium)
- + (constants.m_helion_amu * physics_variables.f_helium3)
+ (constants.M_DEUTERON_AMU * physics_variables.f_deuterium)
+ + (constants.M_TRITON_AMU * physics_variables.f_tritium)
+ + (constants.M_HELION_AMU * physics_variables.f_helium3)
)
# ======================================================================
@@ -3315,16 +3313,16 @@ def plasma_composition() -> None:
# Average atomic masses of injected fuel species in the neutral beams
# Only deuterium and tritium in the beams
physics_variables.m_beam_amu = (
- constants.m_deuteron_amu * (1.0 - current_drive_variables.f_beam_tritium)
- ) + (constants.m_triton_amu * current_drive_variables.f_beam_tritium)
+ constants.M_DEUTERON_AMU * (1.0 - current_drive_variables.f_beam_tritium)
+ ) + (constants.M_TRITON_AMU * current_drive_variables.f_beam_tritium)
# ======================================================================
# Average mass of all ions
physics_variables.m_ions_total_amu = (
(physics_variables.m_fuel_amu * physics_variables.nd_fuel_ions)
- + (constants.m_alpha_amu * physics_variables.nd_alphas)
- + (physics_variables.nd_protons * constants.m_proton_amu)
+ + (constants.M_ALPHA_AMU * physics_variables.nd_alphas)
+ + (physics_variables.nd_protons * constants.M_PROTON_AMU)
+ (physics_variables.m_beam_amu * physics_variables.nd_beam_ions)
)
for imp in range(impurity_radiation_module.N_IMPURITIES):
@@ -3347,30 +3345,30 @@ def plasma_composition() -> None:
(
physics_variables.f_deuterium
* physics_variables.nd_fuel_ions
- / constants.m_deuteron_amu
+ / constants.M_DEUTERON_AMU
)
+ (
physics_variables.f_tritium
* physics_variables.nd_fuel_ions
- / constants.m_triton_amu
+ / constants.M_TRITON_AMU
)
+ (
4.0
* physics_variables.f_helium3
* physics_variables.nd_fuel_ions
- / constants.m_helion_amu
+ / constants.M_HELION_AMU
)
- + (4.0 * physics_variables.nd_alphas / constants.m_alpha_amu)
- + (physics_variables.nd_protons / constants.m_proton_amu)
+ + (4.0 * physics_variables.nd_alphas / constants.M_ALPHA_AMU)
+ + (physics_variables.nd_protons / constants.M_PROTON_AMU)
+ (
(1.0 - current_drive_variables.f_beam_tritium)
* physics_variables.nd_beam_ions
- / constants.m_deuteron_amu
+ / constants.M_DEUTERON_AMU
)
+ (
current_drive_variables.f_beam_tritium
* physics_variables.nd_beam_ions
- / constants.m_triton_amu
+ / constants.M_TRITON_AMU
)
) / physics_variables.dene
for imp in range(impurity_radiation_module.N_IMPURITIES):
@@ -3668,7 +3666,7 @@ def calculate_plasma_current(
# Connor-Hastie asymptotically-correct expression
elif i_plasma_current == 7:
fq = calculate_current_coefficient_hastie(
- alphaj, alphap, bt, triang95, eps, kappa95, p0, constants.rmu0
+ alphaj, alphap, bt, triang95, eps, kappa95, p0, constants.RMU0
)
# Sauter scaling allowing negative triangularity [FED May 2016]
@@ -3687,7 +3685,7 @@ def calculate_plasma_current(
# Main plasma current calculation using the fq value from the different settings
if i_plasma_current != 2:
plasma_current = (
- (constants.twopi / constants.rmu0)
+ (constants.TWOPI / constants.RMU0)
* rminor**2
/ (rmajor * q95)
* fq
@@ -3720,7 +3718,7 @@ def calculate_plasma_current(
kappa,
triang,
len_plasma_poloidal,
- constants.rmu0,
+ constants.RMU0,
)
return bp, qstar, plasma_current
@@ -3782,9 +3780,9 @@ def outplas(self):
# Dimensionless plasma parameters. See reference below.
physics_variables.nu_star = (
1
- / constants.rmu0
- * (15.0e0 * constants.electron_charge**4 * physics_variables.dlamie)
- / (4.0e0 * np.pi**1.5e0 * constants.epsilon0**2)
+ / constants.RMU0
+ * (15.0e0 * constants.ELECTRON_CHARGE**4 * physics_variables.dlamie)
+ / (4.0e0 * np.pi**1.5e0 * constants.EPSILON0**2)
* physics_variables.vol_plasma**2
* physics_variables.rmajor**2
* physics_variables.bt
@@ -3796,7 +3794,7 @@ def outplas(self):
physics_variables.rho_star = np.sqrt(
2.0e0
- * constants.proton_mass
+ * constants.PROTON_MASS
* physics_variables.m_ions_total_amu
* physics_variables.e_plasma_beta
/ (
@@ -3805,7 +3803,7 @@ def outplas(self):
* physics_variables.nd_electron_line
)
) / (
- constants.electron_charge
+ constants.ELECTRON_CHARGE
* physics_variables.bt
* physics_variables.eps
* physics_variables.rmajor
@@ -3814,7 +3812,7 @@ def outplas(self):
physics_variables.beta_mcdonald = (
4.0e0
/ 3.0e0
- * constants.rmu0
+ * constants.RMU0
* physics_variables.e_plasma_beta
/ (physics_variables.vol_plasma * physics_variables.bt**2)
)
@@ -6656,8 +6654,8 @@ def bootstrap_fraction_nevins(
physics_variables.ne0
* physics_variables.te0
* 1.0e3
- * constants.electron_charge
- / (bt**2 / (2.0 * constants.rmu0))
+ * constants.ELECTRON_CHARGE
+ / (bt**2 / (2.0 * constants.RMU0))
)
# Call integration routine using definite integral routine from scipy
@@ -8195,14 +8193,14 @@ def calculate_confinement_time(
# The transport losses is just the electron and ion thermal energies divided by the confinement time.
pden_ion_transport_loss_mw = (
(3 / 2)
- * (constants.electron_charge / 1e3)
+ * (constants.ELECTRON_CHARGE / 1e3)
* nd_ions_total
* tin
/ t_ion_energy_confinement
)
pden_electron_transport_loss_mw = (
(3 / 2)
- * (constants.electron_charge / 1e3)
+ * (constants.ELECTRON_CHARGE / 1e3)
* dene
* ten
/ t_electron_energy_confinement
@@ -8264,17 +8262,17 @@ def calculate_plasma_masses(
"""
# Calculate mass of fuel ions
- m_plasma_fuel_ions = (m_fuel_amu * constants.atomic_mass_unit) * (
+ m_plasma_fuel_ions = (m_fuel_amu * constants.ATOMIC_MASS_UNIT) * (
nd_fuel_ions * vol_plasma
)
- m_plasma_ions_total = (m_ions_total_amu * constants.atomic_mass_unit) * (
+ m_plasma_ions_total = (m_ions_total_amu * constants.ATOMIC_MASS_UNIT) * (
nd_ions_total * vol_plasma
)
- m_plasma_alpha = (nd_alphas * vol_plasma) * constants.alpha_mass
+ m_plasma_alpha = (nd_alphas * vol_plasma) * constants.ALPHA_MASS
- m_plasma_electron = constants.electron_mass * (dene * vol_plasma)
+ m_plasma_electron = constants.ELECTRON_MASS * (dene * vol_plasma)
m_plasma = m_plasma_electron + m_plasma_ions_total
@@ -8312,7 +8310,7 @@ def res_diff_time(rmajor, res_plasma, kappa95):
:param kappa95: plasma elongation at 95% flux surface
"""
- return 2 * constants.rmu0 * rmajor / (res_plasma * kappa95)
+ return 2 * constants.RMU0 * rmajor / (res_plasma * kappa95)
def l_h_threshold_power(
diff --git a/process/physics_functions.py b/process/physics_functions.py
index 4c21d753a8..d61c59dd57 100644
--- a/process/physics_functions.py
+++ b/process/physics_functions.py
@@ -4,8 +4,8 @@
import numpy as np
import process.impurity_radiation as impurity
+from process import constants
from process.data_structure import physics_variables
-from process.fortran import constants
from process.plasma_profiles import PlasmaProfile
logger = logging.getLogger(__name__)
@@ -262,8 +262,8 @@ def fast_alpha_beta(
if physics_variables.f_deuterium < 1.0:
beta_thermal = (
2.0
- * constants.rmu0
- * constants.kiloelectron_volt
+ * constants.RMU0
+ * constants.KILOELECTRON_VOLT
* (dene * ten + nd_ions_total * tin)
/ (bt**2 + bp**2)
)
diff --git a/process/plasma_geometry.py b/process/plasma_geometry.py
index f416d81c68..7cfc39bcb2 100644
--- a/process/plasma_geometry.py
+++ b/process/plasma_geometry.py
@@ -2,15 +2,15 @@
import numpy as np
+from process import constants
from process.data_structure import build_variables, physics_variables
-from process.fortran import constants
logger = logging.getLogger(__name__)
class PlasmaGeom:
def __init__(self):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def plasma_geometry(self) -> None:
"""
@@ -452,7 +452,7 @@ def plasma_volume(
rc = rmajor - rminor + xi
vin = (
- constants.twopi
+ constants.TWOPI
* xi
* (
rc**2 * np.sin(thetai)
@@ -465,7 +465,7 @@ def plasma_volume(
rc = rmajor + rminor - xo
vout = (
- constants.twopi
+ constants.TWOPI
* xo
* (
rc**2 * np.sin(thetao)
diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py
index 4407680a96..cf179e82dc 100644
--- a/process/plasma_profiles.py
+++ b/process/plasma_profiles.py
@@ -4,9 +4,9 @@
import scipy as sp
import process.profiles as profiles
+from process import constants
from process.data_structure import divertor_variables, physics_variables
from process.exceptions import ProcessValueError
-from process.fortran import constants
logger = logging.getLogger(__name__)
@@ -43,7 +43,7 @@ def __init__(self) -> None:
"""
# Default profile_size = 501, but it's possible to experiment with this value.
self.profile_size = 501
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
self.neprofile = profiles.NeProfile(self.profile_size)
self.teprofile = profiles.TeProfile(self.profile_size)
@@ -253,7 +253,7 @@ def calculate_profile_factors(self) -> None:
+ physics_variables.ni0 * physics_variables.ti0
)
* 1.0e3
- * constants.electron_charge
+ * constants.ELECTRON_CHARGE
)
# Pressure profile index (N.B. no pedestal effects included here)
diff --git a/process/power.py b/process/power.py
index a5a3ec6a96..0516735c8a 100644
--- a/process/power.py
+++ b/process/power.py
@@ -3,6 +3,7 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.data_structure import (
build_variables,
@@ -21,14 +22,14 @@
times_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import constants, numerics
+from process.fortran import numerics
logger = logging.getLogger(__name__)
class Power:
def __init__(self):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def pfpwr(self, output: bool):
"""
@@ -752,7 +753,7 @@ def calculate_cryo_loads(self) -> None:
# Calculate electric power requirement for cryogenic plant at tfcoil_variables.temp_tf_cryo (MW)
heat_transport_variables.p_cryo_plant_electric_mw = (
1.0e-6
- * (constants.temp_room - tfcoil_variables.temp_tf_cryo)
+ * (constants.TEMP_ROOM - tfcoil_variables.temp_tf_cryo)
/ (tfcoil_variables.eff_tf_cryo * tfcoil_variables.temp_tf_cryo)
* heat_transport_variables.helpow
)
@@ -774,7 +775,7 @@ def calculate_cryo_loads(self) -> None:
# Calculate electric power requirement for cryogenic plant at tfcoil_variables.tcoolin (MW)
p_tf_cryoal_cryo = (
1.0e-6
- * (constants.temp_room - tfcoil_variables.tcoolin)
+ * (constants.TEMP_ROOM - tfcoil_variables.tcoolin)
/ (tfcoil_variables.eff_tf_cryo * tfcoil_variables.tcoolin)
* heat_transport_variables.helpow_cryal
)
diff --git a/process/process_output.py b/process/process_output.py
index 479778d528..c74645ad49 100644
--- a/process/process_output.py
+++ b/process/process_output.py
@@ -1,26 +1,133 @@
+from contextlib import suppress
+from pathlib import Path
+
import numpy as np
-from process.fortran import constants, process_output_fortran
+from process import constants
+from process.data_structure import global_variables
+from process.fortran import numerics
+
+
+class OutputFileManager:
+ @classmethod
+ def open_files(cls, *, mode="w"):
+ cls._outfile = open( # noqa: SIM115
+ Path(global_variables.output_prefix + "OUT.DAT"), mode
+ )
+ cls._mfile = open( # noqa: SIM115
+ Path(global_variables.output_prefix + "MFILE.DAT"), mode
+ )
+
+ @classmethod
+ def open_idempotence_files(cls):
+ cls._outfile.close()
+ cls._mfile.close()
+
+ cls._outfile = open( # noqa: SIM115
+ Path(global_variables.output_prefix + "IDEM_OUT.DAT"), "w"
+ )
+ cls._mfile = open( # noqa: SIM115
+ Path(global_variables.output_prefix + "IDEM_MFILE.DAT"), "w"
+ )
+
+ @classmethod
+ def close_idempotence_files(cls):
+ Path(cls._outfile.name).unlink()
+ Path(cls._mfile.name).unlink()
+ cls._outfile.close()
+ cls._mfile.close()
+ cls.open_files(mode="a")
+
+ @classmethod
+ def finish(cls):
+ cls._outfile.close()
+ cls._mfile.close()
+
+
+def write(file, string: str):
+ if file == constants.MFILE:
+ OutputFileManager._mfile.write(f"{string}\n") # noqa: SLF001
+ elif file == constants.NOUT:
+ OutputFileManager._outfile.write(f"{string}\n") # noqa: SLF001
+
+
+def ocentr(file, string: str, width: int, *, character="*"):
+ """Write a centred header within a line of characters to a file
+
+ :param file: the integer unit of the file
+ :param string: the heading text
+ :param width: the desired with of the header
+ :param character: the character to pad the heading with (*)
+ """
+ write(file, f"{f' {string} ':{character}^{width}}")
+ write(constants.MFILE, f"# {string} #")
+
+
+def ostars(file, width: int, *, character="*"):
+ """Write a line of characters to a file
+
+ :param file: the integer unit of the file
+ :param width: the desired with of the line
+ :param character: the character to fill the line with (*)
+ """
+ write(file, character * width)
+
+
+def oheadr(file, string: str, *, width: int = 110, character="*"):
+ """Write a centred header within a line of characters between two blank lines
+
+ :param file: the integer unit of the file
+ :param string: the heading text
+ :param width: the desired with of the header
+ :param character: the character to pad the heading with (*)
+ """
+ oblnkl(file)
+ ocentr(file, string, width, character=character)
+ oblnkl(file)
+
+
+def oshead(file, string: str, *, width: int = 80, character="*"):
+ """Write a short centred header within a line of characters between two blank lines
-# necessary to avoid using process_output in the code through
-# two different interfaces
-ocentr = process_output_fortran.ocentr
-ostars = process_output_fortran.ostars
-oheadr = process_output_fortran.oheadr
-oshead = process_output_fortran.oshead
-oblnkl = process_output_fortran.oblnkl
-osubhd = process_output_fortran.osubhd
-ocmmnt = process_output_fortran.ocmmnt
-write = process_output_fortran.write
-dblcol = process_output_fortran.dblcol
-ovarin = process_output_fortran.ovarin
-ovarst = process_output_fortran.ovarst
-obuild = process_output_fortran.obuild
+ :param file: the integer unit of the file
+ :param string: the heading text
+ :param width: the desired with of the header
+ :param character: the character to pad the heading with (*)
+ """
+ oheadr(file, string, width=width, character=character)
+
+
+def oblnkl(file):
+ """Write a blank line to a file
+
+ :param file: the integer unit of the file
+ """
+ write(file, " ")
+
+
+def osubhd(file, string):
+ """Write a subheading between two blank lines
+
+ :param file: the integer unit of the file
+ :param string: the heading text
+ """
+ oblnkl(file)
+ write(file, string)
+ oblnkl(file)
+
+
+def ocmmnt(file, string: str):
+ """Write a comment to a file
+
+ :param file: the integer unit of the file
+ :param string: the comment text
+ """
+ write(file, string)
def ovarre(file, descr: str, varnam: str, value, output_flag: str = ""):
replacement_character = "_"
- if file != constants.mfile:
+ if file != constants.MFILE:
replacement_character = " "
description = f"{descr:<72}".replace(" ", replacement_character)
@@ -28,13 +135,26 @@ def ovarre(file, descr: str, varnam: str, value, output_flag: str = ""):
if isinstance(value, np.ndarray):
value = value.item()
- elif isinstance(value, str):
- value = float(value)
+ if isinstance(value, bytes):
+ # TODO: remove when Fortran is gone
+ value = value.decode().strip()
+ if isinstance(value, str):
+ # try and convert the value to a float
+ # if it fails, leave as a string
+ with suppress(ValueError):
+ value = float(value)
+
+ format_value = f"{value:.17e}" if isinstance(value, float) else f"{value: >12}"
- line = f"{description}{replacement_character} {varname}{replacement_character} {value:.17e} {output_flag}"
+ if varnam.strip("()") in numerics.name_xc:
+ # MDK add ITV label if it is an iteration variable
+ # The ITV flag overwrites the output_flag
+ output_flag = "ITV"
+
+ line = f"{description}{replacement_character} {varname}{replacement_character} {format_value} {output_flag}"
write(file, line)
- if file != constants.mfile:
- ovarre(constants.mfile, descr, varnam, value, output_flag)
+ if file != constants.MFILE:
+ ovarre(constants.MFILE, descr, varnam, value, output_flag)
def ocosts(file, varnam: str, descr: str, value):
@@ -43,3 +163,15 @@ def ocosts(file, varnam: str, descr: str, value):
def ovarrf(file, descr: str, varnam: str, value, output_flag: str = ""):
ovarre(file, descr, varnam, value, output_flag)
+
+
+def ovarin(file, descr: str, varnam: str, value, output_flag: str = ""):
+ ovarre(file, descr, varnam, value, output_flag)
+
+
+def ovarst(file, descr: str, varnam: str, value, output_flag: str = ""):
+ ovarre(file, descr, varnam, value, output_flag)
+
+
+def obuild(file, descr: str, thick: float, total: float, variable_name: str = ""):
+ write(file, f"{descr:<50}{thick:.3e}{' ':<10}{total:.3e} {variable_name}")
diff --git a/process/pulse.py b/process/pulse.py
index 4abc1f956f..00dac052ef 100644
--- a/process/pulse.py
+++ b/process/pulse.py
@@ -1,5 +1,6 @@
import logging
+from process import constants
from process import process_output as po
from process.data_structure import (
constraint_variables,
@@ -9,17 +10,14 @@
pulse_variables,
times_variables,
)
-from process.fortran import (
- constants,
- numerics,
-)
+from process.fortran import numerics
logger = logging.getLogger(__name__)
class Pulse:
def __init__(self):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def run(self, output: bool) -> None:
"""Caller for the pulsed reactor model
diff --git a/process/resistive_tf_coil.py b/process/resistive_tf_coil.py
index fd6eb01f47..a7482a1375 100644
--- a/process/resistive_tf_coil.py
+++ b/process/resistive_tf_coil.py
@@ -3,6 +3,7 @@
import numba
import numpy as np
+from process import constants
from process.data_structure import (
build_variables,
fwbs_variables,
@@ -11,18 +12,16 @@
superconducting_tf_coil_variables,
tfcoil_variables,
)
-from process.fortran import constants
from process.tf_coil import TFCoil
logger = logging.getLogger(__name__)
EPS = np.finfo(1.0).eps
-RMU0 = constants.rmu0
class ResistiveTFCoil(TFCoil):
def __init__(self):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def run(self, output: bool):
"""Run main tfcoil subroutine without outputting."""
diff --git a/process/scan.py b/process/scan.py
index 2cd0a8e3ba..fbeef62675 100644
--- a/process/scan.py
+++ b/process/scan.py
@@ -6,6 +6,7 @@
import process.constraints as constraints
import process.process_output as process_output
+from process import constants
from process.caller import write_output_files
from process.data_structure import (
build_variables,
@@ -15,6 +16,7 @@
current_drive_variables,
divertor_variables,
fwbs_variables,
+ global_variables,
heat_transport_variables,
impurity_radiation_module,
pfcoil_variables,
@@ -24,11 +26,7 @@
tfcoil_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import (
- constants,
- global_variables,
- numerics,
-)
+from process.fortran import numerics
from process.log import logging_model_handler, show_errors
from process.solver_handler import SolverHandler
from process.utilities.f2py_string_patch import (
@@ -178,7 +176,7 @@ def run_scan(self):
# doopt() can also run just an evaluation
ifail = self.doopt()
write_output_files(models=self.models, ifail=ifail)
- show_errors(constants.nout)
+ show_errors(constants.NOUT)
return
if scan_variables.isweep > scan_variables.IPNSCNS:
@@ -207,95 +205,95 @@ def post_optimise(self, ifail: int):
"""
numerics.sqsumsq = (numerics.rcm[: numerics.neqns] ** 2).sum() ** 0.5
- process_output.oheadr(constants.nout, "Numerics")
+ process_output.oheadr(constants.NOUT, "Numerics")
if self.solver == "fsolve":
process_output.ocmmnt(
- constants.nout, "PROCESS has performed an fsolve (evaluation) run."
+ constants.NOUT, "PROCESS has performed an fsolve (evaluation) run."
)
else:
process_output.ocmmnt(
- constants.nout, "PROCESS has performed a VMCON (optimisation) run."
+ constants.NOUT, "PROCESS has performed a VMCON (optimisation) run."
)
if ifail != 1:
- process_output.ovarin(constants.nout, "Error flag", "(ifail)", ifail)
+ process_output.ovarin(constants.NOUT, "Error flag", "(ifail)", ifail)
process_output.oheadr(
- constants.iotty, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION"
+ constants.IOTTY, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION"
)
- process_output.oblnkl(constants.iotty)
+ process_output.oblnkl(constants.IOTTY)
logger.critical(f"Solver returns with ifail /= 1. {ifail=}")
# Error code handler for VMCON
if self.solver == "vmcon":
self.verror(ifail)
- process_output.oblnkl(constants.nout)
- process_output.oblnkl(constants.iotty)
+ process_output.oblnkl(constants.NOUT)
+ process_output.oblnkl(constants.IOTTY)
else:
# Solution found
if self.solver != "fsolve":
process_output.ocmmnt(
- constants.nout, "and found a feasible set of parameters."
+ constants.NOUT, "and found a feasible set of parameters."
)
process_output.oheadr(
- constants.iotty, "PROCESS found a feasible solution"
+ constants.IOTTY, "PROCESS found a feasible solution"
)
else:
process_output.ocmmnt(
- constants.nout, "and found a consistent set of parameters."
+ constants.NOUT, "and found a consistent set of parameters."
)
process_output.oheadr(
- constants.iotty, "PROCESS found a consistent solution"
+ constants.IOTTY, "PROCESS found a consistent solution"
)
- process_output.oblnkl(constants.nout)
- process_output.ovarin(constants.nout, "Error flag", "(ifail)", ifail)
+ process_output.oblnkl(constants.NOUT)
+ process_output.ovarin(constants.NOUT, "Error flag", "(ifail)", ifail)
if numerics.sqsumsq >= 1.0e-2:
- process_output.oblnkl(constants.nout)
+ process_output.oblnkl(constants.NOUT)
process_output.ocmmnt(
- constants.nout,
+ constants.NOUT,
"WARNING: Constraint residues are HIGH; consider re-running",
)
process_output.ocmmnt(
- constants.nout,
+ constants.NOUT,
" with lower values of EPSVMC to confirm convergence...",
)
process_output.ocmmnt(
- constants.nout,
+ constants.NOUT,
" (should be able to get down to about 1.0E-8 okay)",
)
- process_output.oblnkl(constants.nout)
+ process_output.oblnkl(constants.NOUT)
process_output.ocmmnt(
- constants.iotty,
+ constants.IOTTY,
"WARNING: Constraint residues are HIGH; consider re-running",
)
process_output.ocmmnt(
- constants.iotty,
+ constants.IOTTY,
" with lower values of EPSVMC to confirm convergence...",
)
process_output.ocmmnt(
- constants.iotty,
+ constants.IOTTY,
" (should be able to get down to about 1.0E-8 okay)",
)
- process_output.oblnkl(constants.iotty)
+ process_output.oblnkl(constants.IOTTY)
logger.warning(f"High final constraint residues. {numerics.sqsumsq=}")
process_output.ovarin(
- constants.nout, "Number of iteration variables", "(nvar)", numerics.nvar
+ constants.NOUT, "Number of iteration variables", "(nvar)", numerics.nvar
)
process_output.ovarin(
- constants.nout,
+ constants.NOUT,
"Number of constraints (total)",
"(neqns+nineqns)",
numerics.neqns + numerics.nineqns,
)
process_output.ovarin(
- constants.nout, "Optimisation switch", "(ioptimz)", numerics.ioptimz
+ constants.NOUT, "Optimisation switch", "(ioptimz)", numerics.ioptimz
)
# Objective function output: none for fsolve
if self.solver != "fsolve":
process_output.ovarin(
- constants.nout, "Figure of merit switch", "(minmax)", numerics.minmax
+ constants.NOUT, "Figure of merit switch", "(minmax)", numerics.minmax
)
objf_name = string_to_f2py_compatible(
@@ -306,13 +304,13 @@ def post_optimise(self, ifail: int):
numerics.objf_name = objf_name
process_output.ovarst(
- constants.nout,
+ constants.NOUT,
"Objective function name",
"(objf_name)",
numerics.objf_name,
)
process_output.ovarre(
- constants.nout,
+ constants.NOUT,
"Normalised objective function",
"(norm_objf)",
numerics.norm_objf,
@@ -320,7 +318,7 @@ def post_optimise(self, ifail: int):
)
process_output.ovarre(
- constants.nout,
+ constants.NOUT,
"Square root of the sum of squares of the constraint residuals",
"(sqsumsq)",
numerics.sqsumsq,
@@ -328,20 +326,20 @@ def post_optimise(self, ifail: int):
)
if self.solver != "fsolve":
process_output.ovarre(
- constants.nout,
+ constants.NOUT,
"VMCON convergence parameter",
"(convergence_parameter)",
global_variables.convergence_parameter,
"OP ",
)
process_output.ovarin(
- constants.nout,
+ constants.NOUT,
"Number of optimising solver iterations",
"(nviter)",
numerics.nviter,
"OP ",
)
- process_output.oblnkl(constants.nout)
+ process_output.oblnkl(constants.NOUT)
if self.solver == "fsolve":
if ifail == 1:
@@ -349,7 +347,7 @@ def post_optimise(self, ifail: int):
else:
msg = "PROCESS failed to solve using fsolve."
process_output.write(
- constants.nout,
+ constants.NOUT,
f"{msg}\n",
)
else:
@@ -361,7 +359,7 @@ def post_optimise(self, ifail: int):
string2 = "minimise" if numerics.minmax > 0 else "maximise"
process_output.write(
- constants.nout,
+ constants.NOUT,
f"{string1} the optimisation parameters to {string2} the objective function: {objf_name}\n",
)
@@ -383,7 +381,7 @@ def post_optimise(self, ifail: int):
if not written_warning:
written_warning = True
process_output.ocmmnt(
- constants.nout,
+ constants.NOUT,
(
"Certain operating limits have been reached,"
"\n as shown by the following optimisation parameters that are"
@@ -398,14 +396,14 @@ def post_optimise(self, ifail: int):
else:
location, bound = "above", "upper"
process_output.write(
- constants.nout,
+ constants.NOUT,
f" {name:<30}= {xcval} is at or {location} its {bound} bound:"
f" {numerics.itv_scaled_upper_bounds[i] * numerics.scafc[i]}",
)
# Write optimisation parameters to mfile
process_output.ovarre(
- constants.mfile,
+ constants.MFILE,
f2py_compatible_to_string(numerics.lablxc[numerics.ixc[i] - 1]),
f"(itvar{i + 1:03d})",
numerics.xcs[i],
@@ -427,25 +425,25 @@ def post_optimise(self, ifail: int):
)
process_output.ovarre(
- constants.mfile,
+ constants.MFILE,
f"{name} (final value/initial value)",
f"(xcm{i + 1:03d})",
numerics.xcm[i],
)
process_output.ovarre(
- constants.mfile,
+ constants.MFILE,
f"{name} (range normalised)",
f"(nitvar{i + 1:03d})",
xnorm,
)
process_output.ovarre(
- constants.mfile,
+ constants.MFILE,
f"{name} (upper bound)",
f"(boundu{i + 1:03d})",
numerics.itv_scaled_upper_bounds[i] * numerics.scafc[i],
)
process_output.ovarre(
- constants.mfile,
+ constants.MFILE,
f"{name} (lower bound)",
f"(boundl{i + 1:03d})",
numerics.itv_scaled_lower_bounds[i] * numerics.scafc[i],
@@ -453,10 +451,10 @@ def post_optimise(self, ifail: int):
# Write optimisation parameter headings to output file
process_output.osubhd(
- constants.nout, "The solution vector is comprised as follows :"
+ constants.NOUT, "The solution vector is comprised as follows :"
)
process_output.write(
- constants.nout,
+ constants.NOUT,
tabulate(
solution_vector_table,
headers=["", "Final value", "Final / initial"],
@@ -465,7 +463,7 @@ def post_optimise(self, ifail: int):
)
process_output.osubhd(
- constants.nout,
+ constants.NOUT,
"The following equality constraint residues should be close to zero:",
)
@@ -486,7 +484,7 @@ def post_optimise(self, ifail: int):
con1[i],
])
process_output.ovarre(
- constants.mfile,
+ constants.MFILE,
f"{name:<33} normalised residue",
f"(eq_con{numerics.icc[i]:03d})",
con1[i],
@@ -494,7 +492,7 @@ def post_optimise(self, ifail: int):
# Write equality constraints to output file
process_output.write(
- constants.nout,
+ constants.NOUT,
tabulate(
equality_constraint_table,
headers=[
@@ -514,7 +512,7 @@ def post_optimise(self, ifail: int):
# Inequalities not necessarily satisfied when evaluating
if self.solver != "fsolve":
process_output.osubhd(
- constants.nout,
+ constants.NOUT,
"The following inequality constraint residues should be "
"greater than or approximately equal to zero:",
)
@@ -528,14 +526,14 @@ def post_optimise(self, ifail: int):
f"{err[i]} {lab[i]}",
])
process_output.ovarre(
- constants.mfile,
+ constants.MFILE,
f"{name} normalised residue",
f"(ineq_con{numerics.icc[i]:03d})",
numerics.rcm[i],
)
process_output.write(
- constants.nout,
+ constants.NOUT,
tabulate(
inequality_constraint_table,
headers=[
@@ -557,149 +555,149 @@ def verror(self, ifail: int):
an unfeasible result from a VMCON (optimisation) run.
"""
if ifail == -1:
- process_output.ocmmnt(constants.nout, "User-terminated execution of VMCON.")
+ process_output.ocmmnt(constants.NOUT, "User-terminated execution of VMCON.")
process_output.ocmmnt(
- constants.iotty, "User-terminated execution of VMCON."
+ constants.IOTTY, "User-terminated execution of VMCON."
)
elif ifail == 0:
process_output.ocmmnt(
- constants.nout, "Improper input parameters to the VMCON routine."
+ constants.NOUT, "Improper input parameters to the VMCON routine."
)
- process_output.ocmmnt(constants.nout, "PROCESS coding must be checked.")
+ process_output.ocmmnt(constants.NOUT, "PROCESS coding must be checked.")
process_output.ocmmnt(
- constants.iotty, "Improper input parameters to the VMCON routine."
+ constants.IOTTY, "Improper input parameters to the VMCON routine."
)
- process_output.ocmmnt(constants.iotty, "PROCESS coding must be checked.")
+ process_output.ocmmnt(constants.IOTTY, "PROCESS coding must be checked.")
elif ifail == 2:
process_output.ocmmnt(
- constants.nout,
+ constants.NOUT,
"The maximum number of calls has been reached without solution.",
)
process_output.ocmmnt(
- constants.nout,
+ constants.NOUT,
"The code may be stuck in a minimum in the residual space that is significantly above zero.",
)
- process_output.oblnkl(constants.nout)
+ process_output.oblnkl(constants.NOUT)
process_output.ocmmnt(
- constants.nout, "There is either no solution possible, or the code"
+ constants.NOUT, "There is either no solution possible, or the code"
)
process_output.ocmmnt(
- constants.nout, "is failing to escape from a deep local minimum."
+ constants.NOUT, "is failing to escape from a deep local minimum."
)
process_output.ocmmnt(
- constants.nout,
+ constants.NOUT,
"Try changing the variables in IXC, or modify their initial values.",
)
process_output.ocmmnt(
- constants.iotty,
+ constants.IOTTY,
"The maximum number of calls has been reached without solution.",
)
process_output.ocmmnt(
- constants.iotty,
+ constants.IOTTY,
"The code may be stuck in a minimum in the residual space that is significantly above zero.",
)
- process_output.oblnkl(constants.nout)
- process_output.oblnkl(constants.iotty)
+ process_output.oblnkl(constants.NOUT)
+ process_output.oblnkl(constants.IOTTY)
process_output.ocmmnt(
- constants.iotty, "There is either no solution possible, or the code"
+ constants.IOTTY, "There is either no solution possible, or the code"
)
process_output.ocmmnt(
- constants.iotty, "is failing to escape from a deep local minimum."
+ constants.IOTTY, "is failing to escape from a deep local minimum."
)
process_output.ocmmnt(
- constants.iotty,
+ constants.IOTTY,
"Try changing the variables in IXC, or modify their initial values.",
)
elif ifail == 3:
process_output.ocmmnt(
- constants.nout, "The line search required the maximum of 10 calls."
+ constants.NOUT, "The line search required the maximum of 10 calls."
)
process_output.ocmmnt(
- constants.nout, "A feasible solution may be difficult to achieve."
+ constants.NOUT, "A feasible solution may be difficult to achieve."
)
process_output.ocmmnt(
- constants.nout, "Try changing or adding variables to IXC."
+ constants.NOUT, "Try changing or adding variables to IXC."
)
process_output.ocmmnt(
- constants.iotty, "The line search required the maximum of 10 calls."
+ constants.IOTTY, "The line search required the maximum of 10 calls."
)
process_output.ocmmnt(
- constants.iotty, "A feasible solution may be difficult to achieve."
+ constants.IOTTY, "A feasible solution may be difficult to achieve."
)
process_output.ocmmnt(
- constants.iotty, "Try changing or adding variables to IXC."
+ constants.IOTTY, "Try changing or adding variables to IXC."
)
elif ifail == 4:
process_output.ocmmnt(
- constants.nout, "An uphill search direction was found."
+ constants.NOUT, "An uphill search direction was found."
)
process_output.ocmmnt(
- constants.nout, "Try changing the equations in ICC, or"
+ constants.NOUT, "Try changing the equations in ICC, or"
)
- process_output.ocmmnt(constants.nout, "adding new variables to IXC.")
+ process_output.ocmmnt(constants.NOUT, "adding new variables to IXC.")
process_output.ocmmnt(
- constants.iotty, "An uphill search direction was found."
+ constants.IOTTY, "An uphill search direction was found."
)
process_output.ocmmnt(
- constants.iotty, "Try changing the equations in ICC, or"
+ constants.IOTTY, "Try changing the equations in ICC, or"
)
- process_output.ocmmnt(constants.iotty, "adding new variables to IXC.")
+ process_output.ocmmnt(constants.IOTTY, "adding new variables to IXC.")
elif ifail == 5:
process_output.ocmmnt(
- constants.nout, "The quadratic programming technique was unable to"
+ constants.NOUT, "The quadratic programming technique was unable to"
)
- process_output.ocmmnt(constants.nout, "find a feasible point.")
- process_output.oblnkl(constants.nout)
+ process_output.ocmmnt(constants.NOUT, "find a feasible point.")
+ process_output.oblnkl(constants.NOUT)
process_output.ocmmnt(
- constants.nout, "Try changing or adding variables to IXC, or modify"
+ constants.NOUT, "Try changing or adding variables to IXC, or modify"
)
process_output.ocmmnt(
- constants.nout,
+ constants.NOUT,
"their initial values (especially if only 1 optimisation",
)
- process_output.ocmmnt(constants.nout, "iteration was performed).")
+ process_output.ocmmnt(constants.NOUT, "iteration was performed).")
process_output.ocmmnt(
- constants.iotty, "The quadratic programming technique was unable to"
+ constants.IOTTY, "The quadratic programming technique was unable to"
)
- process_output.ocmmnt(constants.iotty, "find a feasible point.")
- process_output.oblnkl(constants.iotty)
+ process_output.ocmmnt(constants.IOTTY, "find a feasible point.")
+ process_output.oblnkl(constants.IOTTY)
process_output.ocmmnt(
- constants.iotty, "Try changing or adding variables to IXC, or modify"
+ constants.IOTTY, "Try changing or adding variables to IXC, or modify"
)
process_output.ocmmnt(
- constants.iotty,
+ constants.IOTTY,
"their initial values (especially if only 1 optimisation",
)
- process_output.ocmmnt(constants.iotty, "iteration was performed).")
+ process_output.ocmmnt(constants.IOTTY, "iteration was performed).")
elif ifail == 6:
process_output.ocmmnt(
- constants.nout, "The quadratic programming technique was restricted"
+ constants.NOUT, "The quadratic programming technique was restricted"
)
process_output.ocmmnt(
- constants.nout, "by an artificial bound, or failed due to a singular"
+ constants.NOUT, "by an artificial bound, or failed due to a singular"
)
- process_output.ocmmnt(constants.nout, "matrix.")
+ process_output.ocmmnt(constants.NOUT, "matrix.")
process_output.ocmmnt(
- constants.nout, "Try changing the equations in ICC, or"
+ constants.NOUT, "Try changing the equations in ICC, or"
)
- process_output.ocmmnt(constants.nout, "adding new variables to IXC.")
+ process_output.ocmmnt(constants.NOUT, "adding new variables to IXC.")
process_output.ocmmnt(
- constants.iotty, "The quadratic programming technique was restricted"
+ constants.IOTTY, "The quadratic programming technique was restricted"
)
process_output.ocmmnt(
- constants.iotty, "by an artificial bound, or failed due to a singular"
+ constants.IOTTY, "by an artificial bound, or failed due to a singular"
)
- process_output.ocmmnt(constants.iotty, "matrix.")
+ process_output.ocmmnt(constants.IOTTY, "matrix.")
process_output.ocmmnt(
- constants.iotty, "Try changing the equations in ICC, or"
+ constants.IOTTY, "Try changing the equations in ICC, or"
)
- process_output.ocmmnt(constants.iotty, "adding new variables to IXC.")
+ process_output.ocmmnt(constants.IOTTY, "adding new variables to IXC.")
def scan_1d(self):
"""Run a 1-D scan."""
@@ -712,7 +710,7 @@ def scan_1d(self):
scan_1d_ifail_dict[iscan] = ifail
write_output_files(models=self.models, ifail=ifail)
- show_errors(constants.nout)
+ show_errors(constants.NOUT)
logging_model_handler.clear_logs()
# outvar now contains results
@@ -767,7 +765,7 @@ def scan_2d(self):
write_output_files(models=self.models, ifail=ifail)
- show_errors(constants.nout)
+ show_errors(constants.NOUT)
logging_model_handler.clear_logs()
scan_2d_ifail_list[iscan_1][iscan_2] = ifail
iscan = iscan + 1
@@ -824,37 +822,37 @@ def scan_2d(self):
def scan_2d_init(self):
process_output.ovarin(
- constants.mfile,
+ constants.MFILE,
"Number of first variable scan points",
"(isweep)",
scan_variables.isweep,
)
process_output.ovarin(
- constants.mfile,
+ constants.MFILE,
"Number of second variable scan points",
"(isweep_2)",
scan_variables.isweep_2,
)
process_output.ovarin(
- constants.mfile,
+ constants.MFILE,
"Scanning first variable number",
"(nsweep)",
scan_variables.nsweep,
)
process_output.ovarin(
- constants.mfile,
+ constants.MFILE,
"Scanning second variable number",
"(nsweep_2)",
scan_variables.nsweep_2,
)
process_output.ovarin(
- constants.mfile,
+ constants.MFILE,
"Scanning second variable number",
"(nsweep_2)",
scan_variables.nsweep_2,
)
process_output.ovarin(
- constants.mfile,
+ constants.MFILE,
"Scanning second variable number",
"(nsweep_2)",
scan_variables.nsweep_2,
@@ -866,22 +864,22 @@ def scan_1d_write_point_header(self, iscan: int):
scan_variables.nsweep, scan_variables.sweep, iscan
)
- process_output.oblnkl(constants.nout)
- process_output.ostars(constants.nout, 110)
+ process_output.oblnkl(constants.NOUT)
+ process_output.ostars(constants.NOUT, 110)
process_output.write(
- constants.nout,
- f"***** Scan point {iscan} of {scan_variables.isweep} : {f2py_compatible_to_string(global_variables.xlabel)}"
- f", {f2py_compatible_to_string(global_variables.vlabel)} = {scan_variables.sweep[iscan - 1]} "
+ constants.NOUT,
+ f"***** Scan point {iscan} of {scan_variables.isweep} : {global_variables.xlabel}"
+ f", {global_variables.vlabel} = {scan_variables.sweep[iscan - 1]} "
"*****",
)
- process_output.ostars(constants.nout, 110)
- process_output.oblnkl(constants.mfile)
- process_output.ovarin(constants.mfile, "Scan point number", "(iscan)", iscan)
+ process_output.ostars(constants.NOUT, 110)
+ process_output.oblnkl(constants.MFILE)
+ process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
print(
f"Starting scan point {iscan} of {scan_variables.isweep} : "
- f"{f2py_compatible_to_string(global_variables.xlabel)} , {f2py_compatible_to_string(global_variables.vlabel)}"
+ f"{global_variables.xlabel} , {global_variables.vlabel}"
f" = {scan_variables.sweep[iscan - 1]}"
)
@@ -898,25 +896,25 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
scan_variables.nsweep_2, scan_variables.sweep_2, iscan_r
)
- process_output.oblnkl(constants.nout)
- process_output.ostars(constants.nout, 110)
+ process_output.oblnkl(constants.NOUT)
+ process_output.ostars(constants.NOUT, 110)
process_output.write(
- constants.nout,
+ constants.NOUT,
f"***** 2D Scan point {iscan} of {scan_variables.isweep * scan_variables.isweep_2} : "
- f"{f2py_compatible_to_string(global_variables.vlabel)} = {scan_variables.sweep[iscan_1 - 1]} and"
- f" {f2py_compatible_to_string(global_variables.vlabel_2)} = {scan_variables.sweep_2[iscan_r - 1]} "
+ f"{global_variables.vlabel} = {scan_variables.sweep[iscan_1 - 1]} and"
+ f" {global_variables.vlabel_2} = {scan_variables.sweep_2[iscan_r - 1]} "
"*****",
)
- process_output.ostars(constants.nout, 110)
- process_output.oblnkl(constants.mfile)
- process_output.ovarin(constants.mfile, "Scan point number", "(iscan)", iscan)
+ process_output.ostars(constants.NOUT, 110)
+ process_output.oblnkl(constants.MFILE)
+ process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
print(
- f"Starting scan point {iscan}: {f2py_compatible_to_string(global_variables.xlabel)}, "
- f"{f2py_compatible_to_string(global_variables.vlabel)} = {scan_variables.sweep[iscan_1 - 1]}"
- f" and {f2py_compatible_to_string(global_variables.xlabel_2)}, "
- f"{f2py_compatible_to_string(global_variables.vlabel_2)} = {scan_variables.sweep_2[iscan_r - 1]} "
+ f"Starting scan point {iscan}: {global_variables.xlabel}, "
+ f"{global_variables.vlabel} = {scan_variables.sweep[iscan_1 - 1]}"
+ f" and {global_variables.xlabel_2}, "
+ f"{global_variables.vlabel_2} = {scan_variables.sweep_2[iscan_r - 1]} "
)
return iscan_r
@@ -924,13 +922,13 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
def scan_1d_write_plot(self):
if scan_variables.first_call_1d:
process_output.ovarin(
- constants.mfile,
+ constants.MFILE,
"Number of scan points",
"(isweep)",
scan_variables.isweep,
)
process_output.ovarin(
- constants.mfile,
+ constants.MFILE,
"Scanning variable number",
"(nsweep)",
scan_variables.nsweep,
diff --git a/process/solver.py b/process/solver.py
index e174d0c3e2..d0efe3c912 100644
--- a/process/solver.py
+++ b/process/solver.py
@@ -15,9 +15,10 @@
)
from scipy.optimize import fsolve
+from process.data_structure import global_variables
from process.evaluators import Evaluators
from process.exceptions import ProcessValueError
-from process.fortran import global_variables, numerics
+from process.fortran import numerics
from process.utilities.f2py_string_patch import f2py_compatible_to_string
logger = logging.getLogger(__name__)
diff --git a/process/stellarator.py b/process/stellarator.py
index 3a0e8984db..cebcefdf6e 100644
--- a/process/stellarator.py
+++ b/process/stellarator.py
@@ -7,9 +7,8 @@
import process.fusion_reactions as reactions
import process.physics_functions as physics_funcs
import process.superconductors as superconductors
-from process import (
- process_output as po,
-)
+from process import constants
+from process import process_output as po
from process.coolprop_interface import FluidProperties
from process.data_structure import (
build_variables,
@@ -18,6 +17,7 @@
current_drive_variables,
divertor_variables,
fwbs_variables,
+ global_variables,
heat_transport_variables,
impurity_radiation_module,
neoclassics_variables,
@@ -32,17 +32,16 @@
times_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import constants, global_variables, numerics
+from process.fortran import numerics
from process.physics import rether
from process.stellarator_config import load_stellarator_config
-from process.utilities.f2py_string_patch import f2py_compatible_to_string
logger = logging.getLogger(__name__)
# NOTE: a different value of electron_charge was used in the original implementation
# making the post-Python results slightly different. As a result, there is a
# relative tolerance on the neoclassics tests of 1e-3
-KEV = 1e3 * constants.electron_charge # Kiloelectron-volt (keV)
+KEV = 1e3 * constants.ELECTRON_CHARGE # Kiloelectron-volt (keV)
class Stellarator:
@@ -90,7 +89,7 @@ def __init__(
:type neoclassics: process.stellarator.Neoclassics
"""
- self.outfile: int = constants.nout
+ self.outfile: int = constants.NOUT
self.first_call_stfwbs = True
self.availability = availability
@@ -192,9 +191,7 @@ def stnewconfig(self):
load_stellarator_config(
stellarator_variables.istell,
- Path(
- f"{f2py_compatible_to_string(global_variables.output_prefix)}stella_conf.json"
- ),
+ Path(f"{global_variables.output_prefix}stella_conf.json"),
)
# If physics_variables.aspect ratio is not in numerics.ixc set it to default value
@@ -859,11 +856,11 @@ def stdiv(self, output: bool):
# Scrape-off temperature in Joules
- e = T_scrape * constants.electron_charge
+ e = T_scrape * constants.ELECTRON_CHARGE
# Sound speed of particles (m/s)
- c_s = np.sqrt(e / (physics_variables.m_fuel_amu * constants.umass))
+ c_s = np.sqrt(e / (physics_variables.m_fuel_amu * constants.UMASS))
# Island size (m)
@@ -4245,8 +4242,8 @@ def stphys(self, output):
physics_variables.beta_fast_alpha
+ physics_variables.beta_beam
+ 2.0e3
- * constants.rmu0
- * constants.electron_charge
+ * constants.RMU0
+ * constants.ELECTRON_CHARGE
* (
physics_variables.dene * physics_variables.ten
+ physics_variables.nd_ions_total * physics_variables.tin
@@ -4258,13 +4255,13 @@ def stphys(self, output):
* physics_variables.beta
* physics_variables.btot
* physics_variables.btot
- / (2.0e0 * constants.rmu0)
+ / (2.0e0 * constants.RMU0)
* physics_variables.vol_plasma
)
physics_variables.rho_star = np.sqrt(
2.0e0
- * constants.proton_mass
+ * constants.PROTON_MASS
* physics_variables.m_ions_total_amu
* physics_variables.e_plasma_beta
/ (
@@ -4273,7 +4270,7 @@ def stphys(self, output):
* physics_variables.nd_electron_line
)
) / (
- constants.electron_charge
+ constants.ELECTRON_CHARGE
* physics_variables.bt
* physics_variables.eps
* physics_variables.rmajor
@@ -4346,14 +4343,14 @@ def stphys(self, output):
physics_variables.fusden_plasma
+ 1.0e6
* physics_variables.p_beam_alpha_mw
- / (constants.dt_alpha_energy)
+ / (constants.DT_ALPHA_ENERGY)
/ physics_variables.vol_plasma
)
physics_variables.fusden_alpha_total = (
physics_variables.fusden_plasma_alpha
+ 1.0e6
* physics_variables.p_beam_alpha_mw
- / (constants.dt_alpha_energy)
+ / (constants.DT_ALPHA_ENERGY)
/ physics_variables.vol_plasma
)
physics_variables.p_dt_total_mw = (
@@ -4935,7 +4932,7 @@ def calc_neoclassics(self):
* impurity_radiation_module.radius_plasma_core_norm
)
dmdt_neo_fuel = (
- dndt_neo_fuel * physics_variables.m_fuel_amu * constants.proton_mass * 1.0e6
+ dndt_neo_fuel * physics_variables.m_fuel_amu * constants.PROTON_MASS * 1.0e6
) # mg
dmdt_neo_fuel_from_e = (
4
@@ -4943,7 +4940,7 @@ def calc_neoclassics(self):
* physics_variables.a_plasma_surface
* impurity_radiation_module.radius_plasma_core_norm
* physics_variables.m_fuel_amu
- * constants.proton_mass
+ * constants.PROTON_MASS
* 1.0e6
) # kg
@@ -5030,7 +5027,7 @@ def st_calc_eff_chi(self):
(
3
* physics_variables.ne0
- * constants.electron_charge
+ * constants.ELECTRON_CHARGE
* physics_variables.te0
* 1e3
* physics_variables.alphat
@@ -5487,12 +5484,12 @@ def neoclassics_calc_KT(self):
def neoclassics_calc_nu(self):
"""Calculates the collision frequency"""
mass = np.array([
- constants.electron_mass,
- constants.proton_mass * 2.0,
- constants.proton_mass * 3.0,
- constants.proton_mass * 4.0,
+ constants.ELECTRON_MASS,
+ constants.PROTON_MASS * 2.0,
+ constants.PROTON_MASS * 3.0,
+ constants.PROTON_MASS * 4.0,
])
- z = np.array([-1.0, 1.0, 1.0, 2.0]) * constants.electron_charge
+ z = np.array([-1.0, 1.0, 1.0, 2.0]) * constants.ELECTRON_CHARGE
# transform the temperature back in eV
# Formula from L. Spitzer.Physics of fully ionized gases. Interscience, New York, 1962
@@ -5501,7 +5498,7 @@ def neoclassics_calc_nu(self):
- 1.15 * np.log10(neoclassics_variables.densities[0])
+ 2.3
* np.log10(
- neoclassics_variables.temperatures[0] / constants.electron_charge
+ neoclassics_variables.temperatures[0] / constants.ELECTRON_CHARGE
)
)
@@ -5544,7 +5541,7 @@ def neoclassics_calc_nu(self):
] + neoclassics_variables.densities[k] * (
z[j] * z[k]
) ** 2 * lnlambda * phixmgx / (
- 4.0 * np.pi * constants.epsilon0**2 * mass[j] ** 2 * v**3
+ 4.0 * np.pi * constants.EPSILON0**2 * mass[j] ** 2 * v**3
)
return neoclassics_calc_nu
@@ -5555,24 +5552,24 @@ def neoclassics_calc_nu_star(self):
kk = (k * neoclassics_variables.temperatures).T
mass = np.array([
- constants.electron_mass,
- constants.proton_mass * 2.0,
- constants.proton_mass * 3.0,
- constants.proton_mass * 4.0,
+ constants.ELECTRON_MASS,
+ constants.PROTON_MASS * 2.0,
+ constants.PROTON_MASS * 3.0,
+ constants.PROTON_MASS * 4.0,
])
v = np.empty((4, self.no_roots))
- v[0, :] = constants.speed_light * np.sqrt(
- 1.0 - (kk[0, :] / (mass[0] * constants.speed_light**2) + 1) ** (-1)
+ v[0, :] = constants.SPEED_LIGHT * np.sqrt(
+ 1.0 - (kk[0, :] / (mass[0] * constants.SPEED_LIGHT**2) + 1) ** (-1)
)
- v[1, :] = constants.speed_light * np.sqrt(
- 1.0 - (kk[1, :] / (mass[1] * constants.speed_light**2) + 1) ** (-1)
+ v[1, :] = constants.SPEED_LIGHT * np.sqrt(
+ 1.0 - (kk[1, :] / (mass[1] * constants.SPEED_LIGHT**2) + 1) ** (-1)
)
- v[2, :] = constants.speed_light * np.sqrt(
- 1.0 - (kk[2, :] / (mass[2] * constants.speed_light**2) + 1) ** (-1)
+ v[2, :] = constants.SPEED_LIGHT * np.sqrt(
+ 1.0 - (kk[2, :] / (mass[2] * constants.SPEED_LIGHT**2) + 1) ** (-1)
)
- v[3, :] = constants.speed_light * np.sqrt(
- 1.0 - (kk[3, :] / (mass[3] * constants.speed_light**2) + 1) ** (-1)
+ v[3, :] = constants.SPEED_LIGHT * np.sqrt(
+ 1.0 - (kk[3, :] / (mass[3] * constants.SPEED_LIGHT**2) + 1) ** (-1)
)
return (
@@ -5600,19 +5597,19 @@ def neoclassics_calc_nu_star_fromT(self, iota):
])
mass = np.array([
- constants.electron_mass,
- constants.proton_mass * 2.0,
- constants.proton_mass * 3.0,
- constants.proton_mass * 4.0,
+ constants.ELECTRON_MASS,
+ constants.PROTON_MASS * 2.0,
+ constants.PROTON_MASS * 3.0,
+ constants.PROTON_MASS * 4.0,
])
- z = np.array([-1.0, 1.0, 1.0, 2.0]) * constants.electron_charge
+ z = np.array([-1.0, 1.0, 1.0, 2.0]) * constants.ELECTRON_CHARGE
# transform the temperature back in eV
# Formula from L. Spitzer.Physics of fully ionized gases. Interscience, New York, 1962
lnlambda = (
32.2
- 1.15 * np.log10(density[0])
- + 2.3 * np.log10(temp[0] / constants.electron_charge)
+ + 2.3 * np.log10(temp[0] / constants.ELECTRON_CHARGE)
)
neoclassics_calc_nu_star_fromT = np.zeros((4,))
@@ -5647,7 +5644,7 @@ def neoclassics_calc_nu_star_fromT(self, iota):
* (z[j] * z[k]) ** 2
* lnlambda
* phixmgx
- / (4.0 * np.pi * constants.epsilon0**2 * mass[j] ** 2 * v**4)
+ / (4.0 * np.pi * constants.EPSILON0**2 * mass[j] ** 2 * v**4)
* physics_variables.rmajor
/ iota
)
@@ -5658,7 +5655,7 @@ def neoclassics_calc_vd(self):
neoclassics_variables.roots
* neoclassics_variables.temperatures[0]
/ (
- constants.electron_charge
+ constants.ELECTRON_CHARGE
* physics_variables.rmajor
* physics_variables.bt
)
@@ -5667,7 +5664,7 @@ def neoclassics_calc_vd(self):
neoclassics_variables.roots
* neoclassics_variables.temperatures[1]
/ (
- constants.electron_charge
+ constants.ELECTRON_CHARGE
* physics_variables.rmajor
* physics_variables.bt
)
@@ -5676,7 +5673,7 @@ def neoclassics_calc_vd(self):
neoclassics_variables.roots
* neoclassics_variables.temperatures[2]
/ (
- constants.electron_charge
+ constants.ELECTRON_CHARGE
* physics_variables.rmajor
* physics_variables.bt
)
@@ -5686,7 +5683,7 @@ def neoclassics_calc_vd(self):
* neoclassics_variables.temperatures[3]
/ (
2.0
- * constants.electron_charge
+ * constants.ELECTRON_CHARGE
* physics_variables.rmajor
* physics_variables.bt
)
@@ -5704,41 +5701,41 @@ def neoclassics_calc_vd(self):
def neoclassics_calc_D11_plateau(self):
"""Calculates the plateau transport coefficients (D11_star sometimes)"""
mass = np.array([
- constants.electron_mass,
- constants.proton_mass * 2.0,
- constants.proton_mass * 3.0,
- constants.proton_mass * 4.0,
+ constants.ELECTRON_MASS,
+ constants.PROTON_MASS * 2.0,
+ constants.PROTON_MASS * 3.0,
+ constants.PROTON_MASS * 4.0,
])
v = np.empty((4, self.no_roots))
- v[0, :] = constants.speed_light * np.sqrt(
+ v[0, :] = constants.SPEED_LIGHT * np.sqrt(
1.0
- (
- neoclassics_variables.kt[0, :] / (mass[0] * constants.speed_light**2)
+ neoclassics_variables.kt[0, :] / (mass[0] * constants.SPEED_LIGHT**2)
+ 1
)
** (-1)
)
- v[1, :] = constants.speed_light * np.sqrt(
+ v[1, :] = constants.SPEED_LIGHT * np.sqrt(
1.0
- (
- neoclassics_variables.kt[1, :] / (mass[1] * constants.speed_light**2)
+ neoclassics_variables.kt[1, :] / (mass[1] * constants.SPEED_LIGHT**2)
+ 1
)
** (-1)
)
- v[2, :] = constants.speed_light * np.sqrt(
+ v[2, :] = constants.SPEED_LIGHT * np.sqrt(
1.0
- (
- neoclassics_variables.kt[2, :] / (mass[2] * constants.speed_light**2)
+ neoclassics_variables.kt[2, :] / (mass[2] * constants.SPEED_LIGHT**2)
+ 1
)
** (-1)
)
- v[3, :] = constants.speed_light * np.sqrt(
+ v[3, :] = constants.SPEED_LIGHT * np.sqrt(
1.0
- (
- neoclassics_variables.kt[3, :] / (mass[3] * constants.speed_light**2)
+ neoclassics_variables.kt[3, :] / (mass[3] * constants.SPEED_LIGHT**2)
+ 1
)
** (-1)
diff --git a/process/structure.py b/process/structure.py
index 0fbf958dca..13069c18f6 100644
--- a/process/structure.py
+++ b/process/structure.py
@@ -3,6 +3,7 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.data_structure import build_variables as bv
from process.data_structure import divertor_variables as divv
@@ -11,7 +12,6 @@
from process.data_structure import physics_variables as pv
from process.data_structure import structure_variables as stv
from process.data_structure import tfcoil_variables as tfv
-from process.fortran import constants
logger = logging.getLogger(__name__)
@@ -26,7 +26,7 @@ class Structure:
"""
def __init__(self) -> None:
- self.outfile = constants.nout # output file unit
+ self.outfile = constants.NOUT # output file unit
def run(self, output: bool = False) -> None:
"""Structure calculation caller
diff --git a/process/superconducting_tf_coil.py b/process/superconducting_tf_coil.py
index 9a56be25b1..69dfef7dc7 100644
--- a/process/superconducting_tf_coil.py
+++ b/process/superconducting_tf_coil.py
@@ -5,12 +5,14 @@
from scipy import optimize
import process.superconductors as superconductors
+from process import constants
from process import process_output as po
from process.data_structure import (
build_variables,
constraint_variables,
divertor_variables,
fwbs_variables,
+ global_variables,
pfcoil_variables,
physics_variables,
rebco_variables,
@@ -18,7 +20,6 @@
tfcoil_variables,
)
from process.exceptions import ProcessValueError
-from process.fortran import constants, global_variables
from process.quench import calculate_quench_protection_current_density
from process.tf_coil import TFCoil
from process.utilities.f2py_string_patch import f2py_compatible_to_string
@@ -38,13 +39,9 @@
}
-RMU0 = constants.rmu0
-EPS = np.finfo(1.0).eps
-
-
class SuperconductingTFCoil(TFCoil):
def __init__(self):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def run(self, output: bool):
"""
@@ -2958,12 +2955,12 @@ def vv_stress_on_quench(
# relevant self-inductances in henry (H)
coil_structure_self_inductance = (
- (constants.rmu0 / np.pi)
+ (constants.RMU0 / np.pi)
* H_coil
* _inductance_factor(H_coil, ri_coil, ro_coil, rm_coil, theta1_coil)
)
vv_self_inductance = (
- (constants.rmu0 / np.pi)
+ (constants.RMU0 / np.pi)
* H_vv
* _inductance_factor(H_vv, ri_vv, ro_vv, rm_vv, theta1_vv)
)
@@ -2990,7 +2987,7 @@ def vv_stress_on_quench(
i2 = (lambda1 / lambda2) * i1
a_vv = (ro_vv + ri_vv) / (ro_vv - ri_vv)
- b_vvi = (constants.rmu0 * (n_tf_coils * n_tf_coil_turns * i0 + i1 + (i2 / 2))) / (
+ b_vvi = (constants.RMU0 * (n_tf_coils * n_tf_coil_turns * i0 + i1 + (i2 / 2))) / (
2 * np.pi * ri_vv
)
j_vvi = i2 / (2 * np.pi * d_vv * ri_vv)
diff --git a/process/tf_coil.py b/process/tf_coil.py
index 8a8ccf5195..8d77fdfb5f 100644
--- a/process/tf_coil.py
+++ b/process/tf_coil.py
@@ -5,12 +5,13 @@
import numba
import numpy as np
-from process import fortran as ft
+from process import constants
from process import process_output as po
from process.build import Build
from process.data_structure import (
build_variables,
fwbs_variables,
+ global_variables,
physics_variables,
rebco_variables,
superconducting_tf_coil_variables,
@@ -18,22 +19,17 @@
)
from process.data_structure import build_variables as bv
from process.exceptions import ProcessValueError
-from process.fortran import constants, global_variables, numerics
-from process.utilities.f2py_string_patch import (
- f2py_compatible_to_string,
-)
+from process.fortran import numerics
logger = logging.getLogger(__name__)
-RMU0 = constants.rmu0
-
class TFCoil:
"""Calculates the parameters of a resistive TF coil system for a fusion power plant"""
def __init__(self, build: Build):
"""Initialise Fortran module variables."""
- self.outfile = ft.constants.nout # output file unit
+ self.outfile = constants.NOUT # output file unit
self.iprint = 0 # switch for writing to output file (1=yes)
self.build = build
self.a_tf_inboard_total = tfcoil_variables.a_tf_inboard_total
@@ -312,7 +308,7 @@ def tf_current(
c_tf_total = (
b_tf_inboard_peak_symmetric
* r_b_tf_inboard_peak
- * (2 * np.pi / constants.rmu0)
+ * (2 * np.pi / constants.RMU0)
)
# Current per TF coil [A]
@@ -665,14 +661,14 @@ def outtf(self):
"OP ",
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
"Inboard leg inner radius (m)",
"(r_tf_inboard_in)",
build_variables.r_tf_inboard_in,
"OP ",
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
"Inboard leg outer radius (m)",
"(r_tf_inboard_out)",
build_variables.r_tf_inboard_out,
@@ -685,28 +681,28 @@ def outtf(self):
tfcoil_variables.i_tf_wp_geom,
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
"Radial position of inner edge and centre of winding pack (m)",
"(r_tf_wp_inboard_inner)",
superconducting_tf_coil_variables.r_tf_wp_inboard_inner,
"OP ",
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
"Radial position of outer edge and of winding pack (m)",
"(r_tf_wp_inboard_outer)",
superconducting_tf_coil_variables.r_tf_wp_inboard_outer,
"OP ",
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
"Radial position of centre of winding pack (m)",
"(r_tf_wp_inboard_centre)",
superconducting_tf_coil_variables.r_tf_wp_inboard_centre,
"OP ",
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
"Minimum toroidal thickness of winding pack (m)",
"(dx_tf_wp_toroidal_min)",
superconducting_tf_coil_variables.dx_tf_wp_toroidal_min,
@@ -820,13 +816,13 @@ def outtf(self):
f" {ii} {tfcoil_variables.r_tf_arc[ii]} {tfcoil_variables.z_tf_arc[ii]}",
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"TF coil arc point {ii} R (m)",
f"(r_tf_arc({ii + 1}))",
tfcoil_variables.r_tf_arc[ii],
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"TF coil arc point {ii} Z (m)",
f"(z_tf_arc({ii + 1}))",
tfcoil_variables.z_tf_arc[ii],
@@ -2175,7 +2171,7 @@ def cntrpst(self):
dcool = 2.0e0 * tfcoil_variables.rcool # Diameter
lcool = 2.0e0 * (bv.z_tf_inside_half + bv.dr_tf_outboard) # Length
tfcoil_variables.ncool = acool / (
- constants.pi * tfcoil_variables.rcool**2
+ constants.PI * tfcoil_variables.rcool**2
) # Number
# Average conductor cross-sectional area to cool (with cooling area)
@@ -2185,7 +2181,7 @@ def cntrpst(self):
/ (bv.z_tf_inside_half + bv.dr_tf_outboard)
+ acool
)
- ro = (acpav / (constants.pi * tfcoil_variables.ncool)) ** 0.5
+ ro = (acpav / (constants.PI * tfcoil_variables.ncool)) ** 0.5
# Inner legs total heating power (to be removed by coolant)
ptot = tfcoil_variables.p_cp_resistive + fwbs_variables.pnuc_cp_tf * 1.0e6
@@ -2198,10 +2194,10 @@ def cntrpst(self):
# --------------
if tfcoil_variables.i_tf_sup == 0:
# Water coolant physical properties
- coolant_density = constants.denh2o
- coolant_cp = constants.cph2o
- coolant_visco = constants.muh2o
- coolant_th_cond = constants.kh2o
+ coolant_density = constants.DENH2O
+ coolant_cp = constants.CPH2O
+ coolant_visco = constants.MUH2O
+ coolant_th_cond = constants.KH2O
# Mass flow rate [kg/s]
cool_mass_flow = acool * coolant_density * tfcoil_variables.vcool
@@ -2274,7 +2270,7 @@ def cntrpst(self):
dtfilmav = ptot / (
h
* 2.0e0
- * constants.pi
+ * constants.PI
* tfcoil_variables.rcool
* tfcoil_variables.ncool
* lcool
@@ -2290,7 +2286,7 @@ def cntrpst(self):
# ******
# Copper conductor
if tfcoil_variables.i_tf_sup == 0:
- conductor_th_cond = constants.k_copper
+ conductor_th_cond = constants.K_COPPER
# Aluminium
elif tfcoil_variables.i_tf_sup == 2:
@@ -3014,7 +3010,7 @@ def tf_coil_self_inductance(
for _ in range(NINTERVALS):
# Field in the dr_bore for unit current
- b = RMU0 / (2.0e0 * np.pi * r)
+ b = constants.RMU0 / (2.0e0 * np.pi * r)
# Find out if there is a dr_bore
if x0 - r < ai:
h_bore = y0 + bi * np.sqrt(1 - ((r - x0) / ai) ** 2)
@@ -3040,7 +3036,7 @@ def tf_coil_self_inductance(
for _ in range(NINTERVALS):
# Field in the dr_bore for unit current
- b = RMU0 / (2.0e0 * np.pi * r)
+ b = constants.RMU0 / (2.0e0 * np.pi * r)
# Find out if there is a dr_bore
if r - x0 < ai:
h_bore = y0 + bi * np.sqrt(1 - ((r - x0) / ai) ** 2)
@@ -3057,8 +3053,8 @@ def tf_coil_self_inductance(
# Picture frame TF coil
ind_tf_coil = (
(z_tf_inside_half + dr_tf_outboard)
- * RMU0
- / constants.pi
+ * constants.RMU0
+ / constants.PI
* np.log(r_tf_outboard_mid / r_tf_inboard_mid)
)
@@ -3088,7 +3084,7 @@ def generic_tf_coil_area_and_masses(self):
tfcoil_variables.tfcryoarea = (
2.0e0
* tfcoil_variables.len_tf_coil
- * constants.twopi
+ * constants.TWOPI
* 0.5e0
* (build_variables.r_tf_inboard_mid + build_variables.r_tf_outboard_mid)
)
@@ -4202,39 +4198,39 @@ def table_format_arrays(a, mult=1, delim="\t\t"):
# MFILE.DAT data
for ii in range(n_tf_bucking + 2):
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"Radial stress at maximum shear of layer {ii + 1} (Pa)",
f"(sig_tf_r_max({ii + 1}))",
sig_tf_r_max[ii],
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"toroidal stress at maximum shear of layer {ii + 1} (Pa)",
f"(sig_tf_t_max({ii + 1}))",
sig_tf_t_max[ii],
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"Vertical stress at maximum shear of layer {ii + 1} (Pa)",
f"(sig_tf_z_max({ii + 1}))",
sig_tf_z_max[ii],
)
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"Von-Mises stress at maximum shear of layer {ii + 1} (Pa)",
f"(sig_tf_vmises_max({ii + 1}))",
sig_tf_vmises_max[ii],
)
if tfcoil_variables.i_tf_tresca == 1 and tfcoil_variables.i_tf_sup == 1:
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"Maximum shear stress for CEA Tresca yield criterion {ii + 1} (Pa)",
f"(s_shear_tf_peak({ii + 1}))",
s_shear_tf_peak[ii],
)
else:
po.ovarre(
- constants.mfile,
+ constants.MFILE,
f"Maximum shear stress for the Tresca yield criterion {ii + 1} (Pa)",
f"(s_shear_tf_peak({ii + 1}))",
s_shear_tf_peak[ii],
@@ -4277,9 +4273,7 @@ def table_format_arrays(a, mult=1, delim="\t\t"):
for k, v in sig_file_data.items()
}
- sig_filename = (
- f2py_compatible_to_string(global_variables.output_prefix) + "SIG_TF.json"
- )
+ sig_filename = global_variables.output_prefix + "SIG_TF.json"
with open(sig_filename, "w") as f:
json.dump(sig_file_data, f)
@@ -4616,10 +4610,10 @@ def extended_plane_strain(
for ii in range(1, nlayers):
currents_enclosed[ii] = currents_enclosed[ii - 1] + currents[ii - 1]
# Factor that multiplies r linearly in the force density
- f_lin_fac[:] = RMU0 / 2.0e0 * d_curr**2
+ f_lin_fac[:] = constants.RMU0 / 2.0e0 * d_curr**2
# Factor that multiplies r reciprocally in the force density
f_rec_fac[:] = (
- RMU0
+ constants.RMU0
/ 2.0e0
* (d_curr * currents_enclosed / np.pi - d_curr**2 * rad[:nlayers] ** 2)
)
@@ -4938,13 +4932,13 @@ def plane_stress(nu, rad, ey, j, nlayers, n_radial_array):
kk = ey / (1 - nu**2)
# Lorentz forces parametrisation coeficients (array equation)
- alpha = 0.5e0 * RMU0 * j**2 / kk
+ alpha = 0.5e0 * constants.RMU0 * j**2 / kk
inner_layer_curr = 0.0e0
for ii in range(nlayers):
beta[ii] = (
0.5e0
- * RMU0
+ * constants.RMU0
* j[ii]
* (inner_layer_curr - np.pi * j[ii] * rad[ii] ** 2)
/ (np.pi * kk[ii])
diff --git a/process/vacuum.py b/process/vacuum.py
index 822473bfdd..b1ac15bfd2 100644
--- a/process/vacuum.py
+++ b/process/vacuum.py
@@ -3,13 +3,13 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.data_structure import build_variables as buv
from process.data_structure import physics_variables as pv
from process.data_structure import tfcoil_variables as tfv
from process.data_structure import times_variables as tv
from process.data_structure import vacuum_variables as vacv
-from process.fortran import constants
logger = logging.getLogger(__name__)
@@ -23,7 +23,7 @@ class Vacuum:
"""
def __init__(self) -> None:
- self.outfile: int = constants.nout
+ self.outfile: int = constants.NOUT
def run(self, output: bool) -> None:
"""Routine to call the vacuum module
@@ -46,7 +46,7 @@ def run(self, output: bool) -> None:
2.0e0
* pv.molflow_plasma_fuelling_required
* pv.m_fuel_amu
- * constants.umass
+ * constants.UMASS
)
self.i_vacuum_pumping = vacv.i_vacuum_pumping
diff --git a/process/water_use.py b/process/water_use.py
index 7c8af2e7ab..bf57fc3bbe 100644
--- a/process/water_use.py
+++ b/process/water_use.py
@@ -1,15 +1,15 @@
import numpy as np
+from process import constants
from process import process_output as po
from process.data_structure import heat_transport_variables, water_usage_variables
-from process.fortran import constants
SECDAY = 86400e0
class WaterUse:
def __init__(self):
- self.outfile = constants.nout
+ self.outfile = constants.NOUT
def run(self, output: bool):
"""
diff --git a/scripts/document_fortran_interface.py b/scripts/document_fortran_interface.py
index 68addb9d77..3ad4f28ff1 100644
--- a/scripts/document_fortran_interface.py
+++ b/scripts/document_fortran_interface.py
@@ -114,7 +114,7 @@ def get_members(
if name[0:2] == "__":
continue
- if type(member) == type(fortran.constants.init_constants): # noqa: E721
+ if type(member) == type(fortran.numerics.init_numerics): # noqa: E721
docstring = member.__doc__
if is_variable(member):
members.append(
diff --git a/source/fortran/constants.f90 b/source/fortran/constants.f90
deleted file mode 100644
index 147b3fb7ba..0000000000
--- a/source/fortran/constants.f90
+++ /dev/null
@@ -1,308 +0,0 @@
-module constants
- !! author: J. Morris (UAKEA)
- !!
- !! Module containing miscellaneous numerical and physical constants
- !!
- !!### References
- !!
- !! -
-#ifndef dp
- use, intrinsic :: iso_fortran_env, only: dp=>real64
-#endif
-
- implicit none
-
- public
-
- ! Standard output unit identifier
- integer, parameter :: iotty = 6
-
- ! Output file unit identifier
- integer, parameter :: nout = 11
-
- ! Plot data file unit identifier
- integer, parameter :: nplot = 12
-
- ! Machine-optimised output file unit
- integer, parameter :: mfile = 13
-
- ! Verbose diagnostics file
- integer, parameter :: vfile = 14
-
- ! Optimisation information output file number
- integer, parameter :: opt_file = 15
-
- ! TF inboard stress radial distributions file number
- integer, parameter :: sig_file = 16
-
- ! degrees to radians, = pi/180
- real(dp), parameter :: degrad = 0.01745329251D0
-
- ! Electron / elementary charge [C]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?e|search_for=electron+charge
- real(dp), parameter :: electron_charge = 1.602176634D-19
-
- ! While the electron charge is a fundamental constant, the electron volt is a derived unit and
- ! is added here for convenience. This allows better syntax and is more readable than using the electron
- ! charge constant directly when working with units of energy.
-
- ! Electron volt [J]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?evj|search_for=electron+volt
- real(dp), parameter :: electron_volt = 1.602176634D-19
-
- ! Kiloelectron volt [J]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?evj|search_for=electron+volt
- real(dp), parameter :: kiloelectron_volt = 1.602176634D-16
-
- ! Unified atomic mass unit [kg]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?u|search_for=atomic+mass+constant
- real(dp), parameter :: atomic_mass_unit = 1.66053906892D-27
-
- ! Electron mass [kg]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?me|search_for=ELECTRON+MASS
- real(dp), parameter :: electron_mass = 9.1093837139D-31
-
- ! Proton mass [kg]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?mp|search_for=PROTON+MASS
- real(dp), parameter :: proton_mass = 1.67262192595D-27
-
- ! Proton mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?mpu
- real(dp), parameter :: m_proton_amu = 1.0072764665789
-
- ! Protium atomic mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=1
- real(dp), parameter :: m_protium_amu = 1.00782503223
-
- ! Deuteron mass [kg]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?md|search_for=DEUTERON+MASS
- real(dp), parameter :: deuteron_mass = 3.3435837768D-27
-
- ! Deuteron mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?mdu
- real(dp), parameter :: m_deuteron_amu = 2.013553212544
-
- ! Triton mass [kg]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?mt|search_for=TRITON+MASS
- real(dp), parameter :: triton_mass = 5.0073567512D-27
-
- ! Triton mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?art
- real(dp), parameter :: m_triton_amu = 3.01550071597
-
- ! Neutron mass [kg]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?mn|search_for=NEUTRON+MASS
- real(dp), parameter :: neutron_mass = 1.67492750056D-27
-
- ! Alpha particle mass [kg]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?mal|search_for=alpha+mass
- real(dp), parameter :: alpha_mass = 6.6446573450D-27
-
- ! Alpha particle mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?malu
- real(dp), parameter :: m_alpha_amu = 4.001506179129
-
- ! Average Helium atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=He
- real(dp), parameter :: m_helium_amu = 4.002602
-
- ! Helion (3He) mass [kg]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?mh|search_for=HELION
- real(dp), parameter :: helion_mass = 5.0064127862D-27
-
- ! Helion (3He) mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?arh
- real(dp), parameter :: m_helion_amu = 3.014932246932
-
- ! Beryllium atom (4Be) mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Be
- real(dp), parameter :: m_beryllium_amu = 9.0121831
-
- ! Average Carbon atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=C
- real(dp), parameter :: m_carbon_amu = 12.0096
-
- ! Average Nitrogen atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=N
- real(dp), parameter :: m_nitrogen_amu = 14.00643
-
- ! Average Oxygen atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=O
- real(dp), parameter :: m_oxygen_amu = 15.99903
-
- ! Average Neon atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Ne
- real(dp), parameter :: m_neon_amu = 20.1797
-
- ! Average Silicon atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Si
- real(dp), parameter :: m_silicon_amu = 28.084
-
- ! Average Argon atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Ar
- real(dp), parameter :: m_argon_amu = 39.948
-
- ! Average Iron atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Fe
- real(dp), parameter :: m_iron_amu = 55.845
-
- ! Average Nickel atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Ni
- real(dp), parameter :: m_nickel_amu = 58.6934
-
- ! Average Krypton atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Kr
- real(dp), parameter :: m_krypton_amu = 83.798
-
- ! Average Xenon atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=Xe
- real(dp), parameter :: m_xenon_amu = 131.293
-
- ! Average Tungsten atom mass [amu]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=W
- real(dp), parameter :: m_tungsten_amu = 183.84
-
- ! Speed of light in vacuum (c) [m/s]
- ! Reference: National Institute of Standards and Technology (NIST)
- ! https://physics.nist.gov/cgi-bin/cuu/Value?c|search_for=light
- real(dp), parameter :: speed_light = 299792458D0
-
- ! Deuterium - Tritium reaction energy [J]
- ! Find the mass difference in the reactancts and products of the D-T reaction
- ! Multiply by the speed of light squared to get the energy released
- real(dp), parameter :: d_t_energy = (((deuteron_mass+triton_mass)-(alpha_mass+neutron_mass))*speed_light**2)
-
- ! Deuterium - Helion (3He) reaction energy [J]
- ! Find the mass difference in the reactancts and products of the D-3He reaction
- ! Multiply by the speed of light squared to get the energy released
- real(dp), parameter :: d_helium_energy = (((deuteron_mass+helion_mass)-(alpha_mass+proton_mass))*speed_light**2)
-
- ! Deuterium - Deuterium (3He producing) reaction energy [J]
- ! Find the mass difference in the reactancts and products of the D-D reaction
- ! Multiply by the speed of light squared to get the energy released
- real(dp), parameter :: dd_helium_energy = (((deuteron_mass+deuteron_mass)-(helion_mass+neutron_mass))*speed_light**2)
-
- ! Deuterium - Deuterium (Triton producing) reaction energy [J]
- ! Find the mass difference in the reactancts and products of the D-D reaction
- ! Multiply by the speed of light squared to get the energy released
- real(dp), parameter :: dd_triton_energy = (((deuteron_mass+deuteron_mass)-(triton_mass+proton_mass))*speed_light**2)
-
- ! Deuterium - Tritium reaction energy fraction carried by neutron
- ! Assuming centre of mass frame as the momenta of the fusion products exceed
- ! those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
- ! Roughly 79.867% of the energy is carried by the neutron
- real(dp), parameter :: dt_neutron_energy_fraction = (alpha_mass/(neutron_mass+alpha_mass))
-
- ! Deuterium - Tritium reaction energy carried by alpha particle neutron [J]
- ! Assuming centre of mass frame as the momenta of the fusion products exceed
- ! those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
- ! Roughly 3.5 MeV of the energy is carried by the alpha particle
- real(dp), parameter :: dt_alpha_energy = (1.0D0-dt_neutron_energy_fraction)*d_t_energy
-
- ! Deuterium - Deuterium (3He producing) reaction energy fraction carried by neutron
- ! Assuming centre of mass frame as the momenta of the fusion products exceed
- ! those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
- ! Roughly 74.935% of the energy is carried by the neutron
- real(dp), parameter :: dd_neutron_energy_fraction = (helion_mass/(neutron_mass+helion_mass))
-
- ! Deuterium - Deuterium (Triton producing) reaction energy fraction carried by proton
- ! Assuming centre of mass frame as the momenta of the fusion products exceed
- ! those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
- ! Roughly 74.960% of the energy is carried by the proton
- real(dp), parameter :: dd_proton_energy_fraction = (triton_mass/(proton_mass+triton_mass))
-
- ! Deuterium - Helion (3He) reaction energy fraction carried by proton
- ! Assuming centre of mass frame as the momenta of the fusion products exceed
- ! those of the fusion reagents by many orders of magnitude. Assumed to be non-relativistic.
- ! Roughly 79.889% of the energy is carried by the proton
- real(dp), parameter :: dhelium_proton_energy_fraction = (alpha_mass/(proton_mass+alpha_mass))
-
- ! Density of Tungsten [kg/m3]
- real(dp), parameter :: den_tungsten = 19250.0D0
-
- ! Room temperature in Kelvin
- ! Assume the room is at 20 degrees Celsius
- real(dp), parameter :: temp_room = 293.15D0
-
- ! pi
- real(dp), parameter :: pi = 3.1415926535897932D0
-
- ! permeability of free space [H/m]
- real(dp), parameter :: rmu0 = 1.256637062D-6
-
- ! 2 pi
- real(dp), parameter :: twopi = 6.2831853071795862D0
-
- ! unified atomic mass unit [kg]
- real(dp), parameter :: umass = 1.660538921D-27
-
- ! permittivity of free space [Farad/m]
- real(dp), parameter :: epsilon0 = 8.85418781D-12
-
- ! specific heat capacity of water (J/kg/K)
- real(dp), parameter :: cph2o = 4180.0D0
-
- ! density of copper (kg/m3)
- real(dp) :: den_copper
-
- ! density of aluminium (kg/m3)
- real(dp) :: den_aluminium
-
- ! density of water (kg/m3)
- real(dp), parameter :: denh2o = 985.0D0
-
- ! Copper thermal conductivity (W/m/K)
- real(dp), parameter :: k_copper = 330.0D0
-
- ! thermal conductivity of water (W/m/K)
- real(dp), parameter :: kh2o = 0.651D0
-
- ! water dynamic viscosity (kg/m/s)
- real(dp), parameter :: muh2o = 4.71D-4
-
- ! Average number of days in a year
- real(dp), parameter :: n_day_year = 365.2425D0
-
- real(dp), parameter :: acceleration_gravity = 9.81D0
- !! Acceleration due to gravity [m/s2]
-
- contains
-
- subroutine init_constants
- !! Initialise module variables
- implicit none
-
- den_copper = 8900.0D0
- den_aluminium = 2700.0D0
- end subroutine init_constants
-end module constants
diff --git a/source/fortran/global_variables.f90 b/source/fortran/global_variables.f90
deleted file mode 100644
index 4d6146aaa7..0000000000
--- a/source/fortran/global_variables.f90
+++ /dev/null
@@ -1,76 +0,0 @@
-module global_variables
- !! author: J. Morris (UKAEA)
- !!
- !! This module contains miscellaneous global variables not well-suited to any
- !! of the other 'variables' modules.
- !!
-#ifndef dp
- use, intrinsic :: iso_fortran_env, only: dp=>real64
-#endif
- implicit none
-
- public
-
- character(len=48) :: icase
- !! power plant type
-
- character(len=180) :: runtitle
- !! short descriptive title for the run
-
- integer :: verbose
- !! switch for turning on/off diagnostic messages
- !!
- !! - =0 turn off diagnostics
- !! - =1 turn on diagnostics
-
- integer :: run_tests
- !! turns on built-in tests if set to 1
-
- integer :: maxcal
- !! maximum number of VMCON iterations
-
- character(len=400) :: fileprefix
- !! input file prefix
-
- character(len=400) :: output_prefix
- !! output file prefix
-
- character(len=25) :: xlabel
- !! scan parameter description label
-
- character(len=25) :: vlabel
- !! scan value name label
-
- character(len=25) :: xlabel_2
- !! scan parameter description label (2nd dimension)
-
- character(len=25) :: vlabel_2
- !! scan value name label (2nd dimension)
-
- integer :: iscan_global
- !! Makes iscan available globally.
-
- real(dp) :: convergence_parameter
- !! VMCON convergence parameter "sum"
-
- contains
-
- subroutine init_global_variables
- !! Initialise global variables
- implicit none
-
- icase = 'Steady-state tokamak model'
- runtitle = "Run Title (change this line using input variable 'runtitle')"
- verbose = 0
- run_tests = 0
- maxcal = 200
- fileprefix = ""
- output_prefix = ""
- xlabel = ""
- vlabel = ""
- xlabel_2 = ""
- vlabel_2 = ""
- iscan_global = 0
- convergence_parameter = 0.0D0
- end subroutine init_global_variables
-end module global_variables
diff --git a/source/fortran/init_module.f90 b/source/fortran/init_module.f90
deleted file mode 100644
index c177a9fe7b..0000000000
--- a/source/fortran/init_module.f90
+++ /dev/null
@@ -1,87 +0,0 @@
-module init_module
-
-#ifndef dp
-use, intrinsic :: iso_fortran_env, only: dp=>real64
-#endif
-
- implicit none
-
- integer, parameter :: nin = 10
-
-contains
-
- subroutine open_files
- use global_variables, only: verbose, fileprefix, output_prefix
- use constants, only: nout, mfile
-
- implicit none
-
- ! Open the input/output external files
- if (trim(fileprefix) == "") then
- open(unit=nin,file="IN.DAT",status='old')
- else
- open(unit=nin,file=trim(fileprefix),status='old')
- end if
-
- open(unit=nout ,file=trim(output_prefix)//'OUT.DAT' ,status='unknown')
- open(unit=mfile ,file=trim(output_prefix)//'MFILE.DAT' ,status='unknown')
-
- end subroutine open_files
-
- subroutine open_idempotence_files
- ! Open new output file and mfile to write output to
- ! This is used when checking model evaluation idempotence, to avoid
- ! polluting the final output file and mfile with intermediate result checks
- use global_variables, only: output_prefix
- use constants, only: nout, mfile
- implicit none
-
- ! Close existing output file and mfile (could be original out and mfiles
- ! or idem scratch files)
- close(unit = nout)
- close(unit = mfile)
- ! Open scratch files with same units
- open(unit=nout, file=trim(output_prefix)//'IDEM_OUT.DAT', action='write', status='replace')
- open(unit=mfile, file=trim(output_prefix)//'IDEM_MFILE.DAT', action='write', status='replace')
- end subroutine open_idempotence_files
-
- subroutine close_idempotence_files
- ! Close the intermediate idempotence-check files, deleting them in the process
- ! Re-open the original OUT.DAT and MFILE.DAT output files, ready to write
- ! the final data, now model evaluation idempotence has been checked
- use global_variables, only: output_prefix
- use constants, only: nout, mfile
- implicit none
-
- ! Close idempotence files, deleting them in the process
- close(unit = nout, status="delete")
- close(unit = mfile, status="delete")
- ! Re-open original output file and mfile, appending future output to them
- open(unit=nout, file=trim(output_prefix)//'OUT.DAT', action='write', position='append')
- open(unit=mfile, file=trim(output_prefix)//'MFILE.DAT', action='write', position='append')
- end subroutine close_idempotence_files
-
- subroutine finish
- ! Originally at the end of the "program", this subroutine writes some final
- ! lines via the output module and then closes any open files. This is
- ! currently called from Python, and will be removed once file handling is
- ! completely dealt with in Python
- ! # TODO Move this output and file handling to Python
-
- use constants, only: iotty, mfile, nout, nplot, opt_file, vfile
- use process_output_fortran, only: oheadr
- use global_variables, only: verbose
- implicit none
-
- call oheadr(nout,'End of PROCESS Output')
- call oheadr(iotty,'End of PROCESS Output')
- call oheadr(nout,'Copy of PROCESS Input Follows')
-
- close(unit = nin)
- close(unit = nout)
- close(unit = nplot)
- close(unit = mfile)
- close(unit = opt_file)
- if (verbose == 1) close(unit = vfile)
- end subroutine finish
-end module init_module
diff --git a/source/fortran/output.f90 b/source/fortran/output.f90
deleted file mode 100755
index 9c7d0f8926..0000000000
--- a/source/fortran/output.f90
+++ /dev/null
@@ -1,748 +0,0 @@
-! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-module process_output_fortran
-
- !! Module containing routines to produce a uniform output style
- !! author: P J Knight, CCFE, Culham Science Centre
- !! N/A
- !! This module contains a number of routines that allow the
- !! program to write output to a file unit in a uniform style.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-#ifndef dp
- use, intrinsic :: iso_fortran_env, only: dp=>real64
-#endif
- implicit none
-
- public
-
-
-
-contains
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine ocentr(file,string,width)
-
- !! Routine to print a centred header within a line of asterisks
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! string : input character string : Character string to be used
- !! width : input integer : Total width of header
- !! This routine writes out a centred header within a line of asterisks.
- !! It cannot cope with a zero-length string; routine
- !! ostars should be used instead.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: active_constraints, ncalls, ipnvars, ioptimz
- use global_variables, only: run_tests, verbose, output_prefix
- use constants, only: mfile
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file, width
- character(len=*), intent(in) :: string
-
- ! Local variables
-
- integer :: lh, nstars, nstars2
- integer, parameter :: maxwidth = 110
- character(len=maxwidth) :: stars
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- stars = repeat('*',maxwidth)
-
- lh = len(string)
-
- if (lh == 0) then
- call ostars(file,width)
- return
- end if
-
- if (width > maxwidth) then
- write(*,*) 'Error in routine OCENTR :'
- write(*,*) 'Maximum width = ',maxwidth
- write(*,*) 'Requested width = ',width
- write(*,*) 'PROCESS stopping.'
- stop 1
- end if
-
- if (lh >= width) then
- write(*,*) 'Error in routine OCENTR :'
- write(*,*) string
- write(*,*) 'This is too long to fit into ',width,' columns.'
- write(*,*) 'PROCESS stopping.'
- stop 1
- end if
-
- ! Number of stars to be printed on the left
-
- nstars = int( (width-lh)/2 ) - 1
-
- ! Number of stars to be printed on the right
-
- nstars2 = width - (nstars+lh+2)
-
- ! Write the whole line
-
- write(file,'(t2,a)') stars(1:nstars)//' '//string//' '//stars(1:nstars2)
-
- write(mfile,'(t2,a)') '#'//' '//string//' '//'#'
-
- end subroutine ocentr
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine ostars(file,width)
-
- !! Routine to print a line of asterisks
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! width : input integer : Total width of header
- !! This routine writes out a line of asterisks.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use global_variables, only: output_prefix, fileprefix
- use constants, only: mfile
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file, width
-
- ! Local variables
-
- integer, parameter :: maxwidth = 110
- character(len=maxwidth) :: stars
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- stars = repeat('*',maxwidth)
-
- write(file,'(1x,a)') stars(1:min(width,maxwidth))
-
- end subroutine ostars
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine oheadr(file,string)
-
- !! Routine to print a centred header within a line of asterisks,
- !! and between two blank lines
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! string : input character string : Character string to be used
- !! This routine writes out a centred header within a line of
- !! asterisks, and between two blank lines.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: active_constraints, sqsumsq, ioptimz
- use global_variables, only: vlabel, run_tests, verbose
- use constants, only: rmu0, pi
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file
- character(len=*), intent(in) :: string
-
- ! Local variables
-
- integer, parameter :: width = 110
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- call oblnkl(file)
- call ocentr(file,string,width)
- call oblnkl(file)
-
- end subroutine oheadr
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine oshead(file,string)
-
- !! Routine to print a short, centred header within a line of asterisks,
- !! and between two blank lines
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! string : input character string : Character string to be used
- !! This routine writes out a short, centred header within a line of
- !! asterisks, and between two blank lines.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use global_variables, only: icase
- use constants, only: pi
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file
- character(len=*), intent(in) :: string
-
- ! Local variables
-
- integer, parameter :: width = 80
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- call oblnkl(file)
- call ocentr(file,string,width)
- call oblnkl(file)
-
- end subroutine oshead
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine oblnkl(file)
-
- !! Routine to print a blank line
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! This routine writes out a simple blank line.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use constants, only: pi, nout
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- write(file,10)
-10 format(' ')
-
- end subroutine oblnkl
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine osubhd(file,string)
-
- !! Routine to print a subheading between two blank lines
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! string : input character string : Character string to be used
- !! This routine writes out a subheading between two blank lines.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use constants, only: iotty, nout
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file
- character(len=*), intent(in) :: string
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- call oblnkl(file)
- call ocmmnt(file,string)
- call oblnkl(file)
-
- end subroutine osubhd
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine ocmmnt(file,string)
-
- !! Routine to print a comment
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! string : input character string : Character string to be used
- !! This routine writes out a comment line.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: boundl, boundu, sqsumsq
- use global_variables, only: icase, vlabel, iscan_global
- use constants, only: rmu0
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file
- character(len=*), intent(in) :: string
-
- ! Local variables
-
- integer, parameter :: maxwidth = 110
- integer :: lh
-! character(len = maxwidth) :: dummy
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- lh = len(trim(string))
-
- if (lh == 0) then
- write(*,*) 'Error in routine OCMMNT :'
- write(*,*) 'A zero-length string is not permitted.'
- write(*,*) 'PROCESS stopping.'
- stop 1
- end if
-
- if (lh >= maxwidth) then
- write(*, *) 'Warning in routine OCMMNT :'
- write(*, '(A)') string
-! write(*,*) 'This is too long to fit into ',maxwidth,' columns.'
- write(*, '(A,i3,A)') 'This is longer than ',maxwidth,' columns.' ! MK 28/10/2016 Modified previous output to reflect warning message
- !write(*,*) 'PROCESS stopping.'
- !stop 1
- end if
-! dummy = trim(string)
- write(file,'(t2,a)') trim(string)
- !write(file,'(t2,a)') dummy
- end subroutine ocmmnt
-
- subroutine write(file, string)
- !! Write a string to file.
- !! file : input integer : Fortran output unit identifier
- !! string : input character string : Character string to be used
- implicit none
-
- ! Arguments
- integer, intent(in) :: file
- character(len=*), intent(in) :: string
-
- write(file,*) trim(string)
- end subroutine write
-
- subroutine dblcol(file, desc, val1, val2)
- !! Write a description and 2 columns of values to 2dp in standard notation.
- !! file : input integer : Fortran output unit identifier
- !! desc : input character string : Character string to be used
- !! val1 : input real : Value of the left variable
- !! val2 : input real : Value of the right variable
- implicit none
-
- ! Arguments
- integer, intent(in) :: file
- character(len=70), intent(in) :: desc
- real(8), intent(in) :: val1, val2
-
- write(file,10) desc, val1, val2
- 10 format(1x,a,t75,f10.2,t100,f10.2)
- end subroutine dblcol
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine ovarre(file,descr,varnam,value,output_flag)
-
- !! Routine to print out the details of a floating-point
- !! variable using 'E' format
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! descr : input character string : Description of the variable
- !! varnam : input character string : Name of the variable
- !! value : input real : Value of the variable
- !! output_flag : optional character
- !! This routine writes out the description, name and value of a
- !! double precision variable in E format (e.g.
- !! -1.234E+04).
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: name_xc
- use global_variables, only: icase, vlabel
- use constants, only: mfile, nout
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file
- character(len=*), intent(in) :: descr, varnam
- real(8), intent(in) :: value
- character(len=3), intent(in), optional :: output_flag
-
- ! Local variables
-
- character(len=72) :: dum72
- character(len=30) :: dum20
- character(len=20) :: stripped
- character(len=3) :: flag
- integer :: dotindex
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- ! Replace descr and varnam with dummy strings of the correct length.
- ! This counters problems that would occur if the two original strings
- ! were the wrong length.
-
- dum72 = descr
- dum20 = varnam
- ! Remove the "(" and ")" from the varnam
- stripped = varnam(2:len(varnam)-1)
-
- ! May need to strip Python module name (e.g. the pfv. from pfv.j_cs_flat_top_end)
- ! This ensures the ITV flag is still added when required in output files
- dotindex = scan(stripped,".")
- stripped = stripped(dotindex+1:)
-
- if (present(output_flag)) then
- flag = output_flag
- else
- flag = ''
- end if
-
- if (any(name_xc == stripped)) flag = 'ITV'
-
- if (file /= mfile) then
- ! MDK add ITV label if it is an iteration variable
- ! The ITV flag overwrites the output_flag
- write(file,20) dum72, dum20, value, flag
- end if
-
- call underscore(dum72)
- call underscore(dum20)
- write(mfile,10) dum72, dum20, value, flag
-
- ! MFILE.DAT format
- ! Machine epsilon for double ~2.22e-16, hence require 17 sig figs in significand
- ! for full precision of a double float
-10 format(1x,a,t75,a,t110,ES23.16e2," ",a,t10)
- ! OUT.DAT format
-20 format(1x,a,t75,a,t100,1pe10.3, t112, a)
-
- end subroutine ovarre
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine ovarin(file,descr,varnam,value,output_flag)
-
- !! Routine to print out the details of an integer variable
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! descr : input character string : Description of the variable
- !! varnam : input character string : Name of the variable
- !! value : input integer : Value of the variable
- !! output_flag : optional character
- !! This routine writes out the description, name and value of an
- !! integer variable.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: name_xc, icc, ioptimz
- use global_variables, only: xlabel_2, iscan_global
- use constants, only: mfile, nout
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file
- character(len=*), intent(in) :: descr, varnam
- integer, intent(in) :: value
- character(len=3), intent(in), optional :: output_flag
-
- ! Local variables
-
- character(len=72) :: dum72
- character(len=30) :: dum20
- character(len=20) :: stripped
- character(len=3) :: flag
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- ! Replace descr and varnam with dummy strings of the correct length.
- ! This counters problems that would occur if the two original strings
- ! were the wrong length.
-
- dum72 = descr
- dum20 = varnam
- stripped = varnam(2:len(varnam)-1)
- if (present(output_flag)) then
- flag = output_flag
- else
- flag = ''
- end if
-
- if (any(name_xc == stripped)) flag = 'ITV'
-
- if (file /= mfile) then
- ! MDK add ITV label if it is an iteration variable
- ! The ITV flag overwrites the output_flag
- write(file,20) dum72, dum20, value, flag
- end if
-
- call underscore(dum72)
- call underscore(dum20)
- write(mfile,10) dum72, dum20, value, flag
-
-10 format(1x,a,t75,a,t110,i10," ",a,t10)
-20 format(1x,a,t75,a,t100,i10,t112, a)
-
- end subroutine ovarin
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine ovarst(file,descr,varnam,value)
-
- !! Routine to print out the details of a character variable
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! descr : input character string : Description of the variable
- !! varnam : input character string : Name of the variable
- !! value : input character string : Value of the variable
- !! This routine writes out the description, name and value of a
- !! character string variable.
- !! None
- !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: sqsumsq
- use constants, only: mfile
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file
- character(len=*), intent(in) :: descr, varnam, value
-
- ! Local variables
-
- character(len=72) :: dum72
- character(len=30) :: dum20
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- ! Replace descr and varnam with dummy strings of the correct length.
- ! This counters problems that would occur if the two original strings
- ! were the wrong length.
-
- dum72 = descr
- dum20 = varnam
-
- if (file /= mfile) then
- write(file,10) dum72, dum20, value
- end if
-
- call underscore(dum72)
- call underscore(dum20)
- write(mfile,10) dum72, dum20, value
-
-10 format(1x,a,t75,a,t110,a)
-
- end subroutine ovarst
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine obuild(file,descr,thick,total,variable_name)
-
- !! Routine to print out a description, the thickness and
- !! summed build of a component of the radial or vertical build
- !! author: P J Knight, CCFE, Culham Science Centre
- !! file : input integer : Fortran output unit identifier
- !! descr : input character string : Description of the component
- !! thick : input real : Thickness of the component (m)
- !! total : input real : Total build, including this component (m)
- !! This routine writes out a description, the thickness and
- !! summed build of a component of the radial or vertical build.
- !! !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: boundl, boundu
- use constants, only: electron_charge
- implicit none
-
- ! Arguments
-
- integer, intent(in) :: file
- character(len=*), intent(in) :: descr
- character(len=*), optional :: variable_name
- real(8), intent(in) :: thick, total
-
- ! Local variables
-
- character(len=50) :: dum30
- character(len=40) :: dum20
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- ! Replace descr with dummy string of the correct length.
- ! This counters problems that would occur if the original string
- ! was the wrong length.
-
- dum30 = descr
- if (present(variable_name)) then
- dum20 = variable_name
- else
- dum20 = ''
- end if
-
- write(file,10) dum30, thick, total, dum20
-10 format(1x,a,t42,f10.3,t58,f10.3,t71,a)
-
- end subroutine obuild
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- subroutine underscore(string)
-
- !! Routine that converts spaces in a string to underscores
- !! author: P J Knight, CCFE, Culham Science Centre
- !! string : input/output string : character string of interest
- !! This routine converts any space characters in the string
- !! to underscore characters.
- !! None
- !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: active_constraints, boundu, boundl
- use constants, only: electron_charge
- implicit none
-
- ! Arguments
-
- character(len=*), intent(inout) :: string
-
- ! Local variables
-
- integer :: loop, i
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- i = index(string, ' ')
- if (i > 0) then
- do loop = i, len(string)
- if (string(loop:loop) == ' ') string(loop:loop) = '_'
- end do
- end if
-
- end subroutine underscore
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- function int2char(i)
-
- !! Converts a single-digit integer into a character string
- !! author: P J Knight, CCFE, Culham Science Centre
- !! i : input integer : must be between 0 and 9
- !! This is a very simple routine that converts a single-digit
- !! integer into a character string. If the integer is outside
- !! the range 0 to 9 the program stops with an error.
- !! None
- !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: epsvmc, boundu
- use constants, only: rmu0
- implicit none
-
- character(len=1) :: int2char
-
- ! Arguments
-
- integer, intent(in) :: i
-
- ! Local variables
-
- character(len=10), parameter :: number = '0123456789'
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- if ((i < 0).or.(i > 9)) then
- write(*,*) 'INT2CHAR: illegal argument'
- stop 1
- end if
-
- int2char = number(i+1:i+1)
-
- end function int2char
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- function int_to_string2(i)
-
- !! Converts a positive integer into a two-digit character string
- !! author: P J Knight, CCFE, Culham Science Centre
- !! i : input integer : must be between 0 and 99
- !! This routine converts a positive integer into a two-digit
- !! character string.
- !! If the integer is negative, the routine stops with an error.
- !! If the integer is greater than 99, the routine returns a
- !! string containing its last two digits.
- !! None
- !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: boundu
- use constants, only: pi
- implicit none
-
- character(len=2) :: int_to_string2
-
- ! Arguments
-
- integer, intent(in) :: i
-
- ! Local variables
-
- character(len=1) :: a0, a1
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- if (i < 0) then
- write(*,*) 'INT_TO_STRING2: illegal argument'
- stop 1
- end if
-
- a0 = int2char(mod(i,10))
- a1 = int2char(mod(int(i/10),10))
-
- int_to_string2 = a1//a0
-
- end function int_to_string2
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- function int_to_string3(i)
-
- !! Converts a positive integer into a 3-digit character string
- !! author: P J Knight, CCFE, Culham Science Centre
- !! i : input integer : must be between 0 and 99
- !! This routine converts a positive integer into a three-digit
- !! character string.
- !! If the integer is negative, the routine stops with an error.
- !! If the integer is greater than 999, the routine returns a
- !! string containing its last three digits.
- !! None
- !
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use numerics, only: boundu
- use constants, only: pi
- implicit none
-
- character(len=3) :: int_to_string3
-
- ! Arguments
-
- integer, intent(in) :: i
-
- ! Local variables
-
- character(len=1) :: a0, a1, a2
-
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- if (i < 0) then
- write(*,*) 'INT_TO_STRING3: illegal argument'
- stop 1
- end if
-
- a0 = int2char(mod(i,10))
- a1 = int2char(mod(int(i/10),10))
- a2 = int2char(mod(int(i/100),100))
-
- int_to_string3 = a2//a1//a0
-
- end function int_to_string3
-
-end module process_output_fortran
diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py
index 9d82a7c17e..bb2c799bcc 100644
--- a/tests/integration/test_pfcoil_int.py
+++ b/tests/integration/test_pfcoil_int.py
@@ -13,6 +13,7 @@
import pytest
from numpy.testing import assert_array_almost_equal
+from process import constants
from process.cs_fatigue import CsFatigue
from process.data_structure import build_variables as bv
from process.data_structure import fwbs_variables as fwbsv
@@ -20,7 +21,6 @@
from process.data_structure import physics_variables as pv
from process.data_structure import tfcoil_variables as tfv
from process.data_structure import times_variables as tv
-from process.fortran import constants
from process.init import init_all_module_vars
from process.pfcoil import PFCoil, fixb, mtrx
diff --git a/tests/unit/test_ccfe_hcpb.py b/tests/unit/test_ccfe_hcpb.py
index 46ee44306e..d45e467467 100644
--- a/tests/unit/test_ccfe_hcpb.py
+++ b/tests/unit/test_ccfe_hcpb.py
@@ -8,14 +8,12 @@
current_drive_variables,
divertor_variables,
fwbs_variables,
+ global_variables,
heat_transport_variables,
physics_variables,
primary_pumping_variables,
tfcoil_variables,
)
-from process.fortran import (
- global_variables,
-)
from process.fw import Fw
from process.hcpb import CCFE_HCPB
diff --git a/tests/unit/test_input.py b/tests/unit/test_input.py
index 47dcccf7bb..ed36838dd3 100644
--- a/tests/unit/test_input.py
+++ b/tests/unit/test_input.py
@@ -1,14 +1,11 @@
import numpy as np
import pytest
+import process.data_structure as data_structure
import process.init as init
import process.input as process_input
from process import fortran
from process.exceptions import ProcessValidationError
-from process.utilities.f2py_string_patch import (
- f2py_compatible_to_string,
- string_to_f2py_compatible,
-)
def _create_input_file(directory, content: str):
@@ -62,9 +59,8 @@ def test_parse_real(epsvmc, expected, tmp_path):
Program to get the expected value for 0.008 provided at https://github.com/ukaea/PROCESS/pull/3067
"""
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- _create_input_file(tmp_path, f"epsvmc = {epsvmc}"),
+ data_structure.global_variables.fileprefix = _create_input_file(
+ tmp_path, f"epsvmc = {epsvmc}"
)
init.init_process()
@@ -88,9 +84,8 @@ def test_exact_parsing(value, tmp_path):
These tests failed using the old input parser and serve to show that the Python parser generally
produces more accurate floats and accumulates less error.
"""
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- _create_input_file(tmp_path, f"epsvmc = {value}"),
+ data_structure.global_variables.fileprefix = _create_input_file(
+ tmp_path, f"epsvmc = {value}"
)
init.init_process()
@@ -98,27 +93,21 @@ def test_exact_parsing(value, tmp_path):
def test_parse_input(tmp_path):
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- _create_input_file(
- tmp_path,
- ("runtitle = my run title\nioptimz = -2\nepsvmc = 0.6\nboundl(1) = 0.5"),
- ),
+ data_structure.global_variables.fileprefix = _create_input_file(
+ tmp_path,
+ ("runtitle = my run title\nioptimz = -2\nepsvmc = 0.6\nboundl(1) = 0.5"),
)
init.init_process()
- assert (
- f2py_compatible_to_string(fortran.global_variables.runtitle) == "my run title"
- )
+ assert data_structure.global_variables.runtitle == "my run title"
assert fortran.numerics.ioptimz == -2
assert pytest.approx(fortran.numerics.epsvmc) == 0.6
assert pytest.approx(fortran.numerics.boundl[0]) == 0.5
def test_input_choices(tmp_path):
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- _create_input_file(tmp_path, ("ioptimz = -1")),
+ data_structure.global_variables.fileprefix = _create_input_file(
+ tmp_path, ("ioptimz = -1")
)
with pytest.raises(ProcessValidationError):
@@ -129,9 +118,8 @@ def test_input_choices(tmp_path):
("input_file_value"), ((-0.01,), (1.1,)), ids=("violate lower", "violate upper")
)
def test_input_range(tmp_path, input_file_value):
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- _create_input_file(tmp_path, (f"epsvmc = {input_file_value}")),
+ data_structure.global_variables.fileprefix = _create_input_file(
+ tmp_path, (f"epsvmc = {input_file_value}")
)
# check that the test data doesn't change
@@ -142,9 +130,8 @@ def test_input_range(tmp_path, input_file_value):
def test_input_array_when_not(tmp_path):
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- _create_input_file(tmp_path, ("epsvmc(1) = 0.5")),
+ data_structure.global_variables.fileprefix = _create_input_file(
+ tmp_path, ("epsvmc(1) = 0.5")
)
with pytest.raises(ProcessValidationError):
@@ -152,9 +139,8 @@ def test_input_array_when_not(tmp_path):
def test_input_not_array_when_is(tmp_path):
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- _create_input_file(tmp_path, ("boundl = 0.5")),
+ data_structure.global_variables.fileprefix = _create_input_file(
+ tmp_path, ("boundl = 0.5")
)
with pytest.raises(ProcessValidationError):
@@ -162,9 +148,8 @@ def test_input_not_array_when_is(tmp_path):
def test_input_float_when_int(tmp_path):
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- _create_input_file(tmp_path, ("ioptimz = 0.5")),
+ data_structure.global_variables.fileprefix = _create_input_file(
+ tmp_path, ("ioptimz = 0.5")
)
with pytest.raises(ProcessValidationError):
@@ -172,10 +157,10 @@ def test_input_float_when_int(tmp_path):
def test_input_array(tmp_path):
- fortran.global_variables.fileprefix = string_to_f2py_compatible(
- fortran.global_variables.fileprefix,
- _create_input_file(tmp_path, ("boundl = 0.1, 0.2, 1.0, 0.0, 1.0e2")),
+ data_structure.global_variables.fileprefix = _create_input_file(
+ tmp_path, ("boundl = 0.1, 0.2, 1.0, 0.0, 1.0e2")
)
+
init.init_process()
np.testing.assert_array_equal(
fortran.numerics.boundl[:6], [0.1, 0.2, 1.0, 0.0, 1.0e2, 0]
diff --git a/tests/unit/test_main.py b/tests/unit/test_main.py
index 98785f1838..9360af8f7e 100644
--- a/tests/unit/test_main.py
+++ b/tests/unit/test_main.py
@@ -6,11 +6,8 @@
import pytest
-from process import fortran, main
+from process import data_structure, main
from process.main import Process, SingleRun, VaryRun
-from process.utilities.f2py_string_patch import (
- f2py_compatible_to_string,
-)
def test_main(monkeypatch):
@@ -136,6 +133,7 @@ def single_run(monkeypatch, input_file, tmp_path):
single_run = SingleRun()
temp_input_file = shutil.copy(input_file, tmp_path / Path(input_file).name)
+ print(temp_input_file)
single_run.input_file = str(temp_input_file)
single_run.models = None
@@ -168,44 +166,36 @@ def test_set_input(single_run, monkeypatch, input_file):
monkeypatch.setattr(single_run, "input_file", input_file, raising=False)
# Mocking undo trys to set the value as none
- # TODO Create our own monkeypatch for strings (if needed)
- # Mock the Fortran set
- # monkeypatch.setattr(fortran.global_variables, "fileprefix", string_to_f2py_compatible(fortran.global_variables.fileprefix,None))
- # fortran.global_variables.test_setting_string()
# Mocks set up, can now run set_input()
single_run.set_input()
- # Check path has been set in the Fortran (mocked above)
- result = f2py_compatible_to_string(fortran.global_variables.fileprefix)
- assert result == expected
+ # Check path has been set
+ assert data_structure.global_variables.fileprefix == expected
def test_set_output(single_run, monkeypatch):
- """Check output filename setting in the Fortran.
+ """Check output filename set correctly.
:param single_run: single_run fixture
:type single_run: SingleRun
:param monkeypatch: monkeypatch fixture
:type monkeypatch: object
"""
- # Expected output prefix stored in Fortran
+ # Expected output prefix
expected = "output_prefix"
# Mock self.filename_prefix on single_run with the value of expected
monkeypatch.setattr(single_run, "filename_prefix", expected, raising=False)
# Mocking undo trys to set the value as none
- # TODO Create our own monkeypatch for strings (if needed)
- # Mock the Fortran set
- # monkeypatch.setattr(fortran.global_variables, "output_prefix", None)
- # Run the method, and extract the value from the Fortran
+ # monkeypatch.setattr(data_structure.global_variables, "output_prefix", None)
+ # Run the method, and extract the value
single_run.set_output()
- # Convert string from byte-string for comparison
- result = f2py_compatible_to_string(fortran.global_variables.output_prefix)
- assert result == expected
+
+ assert data_structure.global_variables.output_prefix == expected
def test_initialise(single_run, monkeypatch):
- """Test that the init_module can be called in the Fortran.
+ """Test that the init_module runs without crashing
:param single_run: single_run fixture
:type single_run: SingleRun
diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py
index eb5e4f10e2..4129e1479a 100644
--- a/tests/unit/test_physics.py
+++ b/tests/unit/test_physics.py
@@ -5,6 +5,7 @@
import numpy as np
import pytest
+from process import constants
from process.current_drive import (
CurrentDrive,
ElectronBernstein,
@@ -18,7 +19,6 @@
impurity_radiation_module,
physics_variables,
)
-from process.fortran import constants
from process.impurity_radiation import initialise_imprad
from process.physics import (
Physics,
@@ -1241,7 +1241,7 @@ def test_calculate_plasma_current_peng(arguments, expected):
"kappa": 1.85,
"delta": 0.5,
"perim": 24,
- "rmu0": constants.rmu0,
+ "rmu0": constants.RMU0,
},
3.4401302177092803,
),
@@ -1256,7 +1256,7 @@ def test_calculate_plasma_current_peng(arguments, expected):
"kappa": 1.85,
"delta": 0.5,
"perim": 24,
- "rmu0": constants.rmu0,
+ "rmu0": constants.RMU0,
},
2.9310284627233205,
),
@@ -1271,7 +1271,7 @@ def test_calculate_plasma_current_peng(arguments, expected):
"kappa": 1.85,
"delta": 0.5,
"perim": 24,
- "rmu0": constants.rmu0,
+ "rmu0": constants.RMU0,
},
0.8377580413333333,
),
@@ -1287,7 +1287,7 @@ def test_calculate_beta_limit():
def test_conhas():
assert calculate_current_coefficient_hastie(
- 5, 5, 12, 0.5, 0.33, 1.85, 2e3, constants.rmu0
+ 5, 5, 12, 0.5, 0.33, 1.85, 2e3, constants.RMU0
) == pytest.approx(2.518876726889116)
diff --git a/tests/unit/test_sctfcoil.py b/tests/unit/test_sctfcoil.py
index 2bf9552bb8..09029cf192 100644
--- a/tests/unit/test_sctfcoil.py
+++ b/tests/unit/test_sctfcoil.py
@@ -9,13 +9,11 @@
constraint_variables,
divertor_variables,
fwbs_variables,
+ global_variables,
physics_variables,
superconducting_tf_coil_variables,
tfcoil_variables,
)
-from process.fortran import (
- global_variables,
-)
from process.superconducting_tf_coil import SuperconductingTFCoil
diff --git a/tracking/tracking_data.py b/tracking/tracking_data.py
index 59208ca048..97a2aa6165 100644
--- a/tracking/tracking_data.py
+++ b/tracking/tracking_data.py
@@ -559,10 +559,10 @@ def _get_variables(cls, fortran_module):
variables = []
for name, function in inspect.getmembers(fortran_module):
- # type(data_structure.physics_variables.constants.init_constants) => fortran subroutine
+ # type(data_structure.physics_variables.numerics.init_numerics) => fortran subroutine
# if its not a fortran subroutine, its a variable
# because type `fortran` cannot be checked as a type
- if type(function) != type(fortran.constants.init_constants): # noqa: E721
+ if type(function) != type(fortran.numerics.init_numerics): # noqa: E721
variables.append(name)
return variables