Fix Stellarator to allow it to run on open-source process

process/main.py

@@ -347,6 +347,7 @@ def __init__(self, input_file, solver="vmcon"):
         self.validate_input()
         self.init_module_vars()
         self.set_filenames()
+        self.initialise()
         self.models = Models()
         self.solver = solver
 
@@ -355,8 +356,6 @@ def run(self):
 
         This is separate from init to allow model instances to be modified before a run.
         """
-        self.set_filenames()
-        self.initialise()
         self.validate_user_model()
         self.run_tests()
         self.call_solver()
@@ -603,6 +602,9 @@ def __init__(self):
         self.blanket_library = BlanketLibrary(fw=self.fw)
         self.ccfe_hcpb = CCFE_HCPB(blanket_library=self.blanket_library)
         self.current_drive = CurrentDrive()
+        self.physics = Physics(
+            plasma_profile=self.plasma_profile, current_drive=self.current_drive
+        )
         self.stellarator = Stellarator(
             availability=self.availability,
             buildings=self.buildings,
@@ -612,9 +614,7 @@ def __init__(self):
             plasma_profile=self.plasma_profile,
             hcpb=self.ccfe_hcpb,
             current_drive=self.current_drive,
-        )
-        self.physics = Physics(
-            plasma_profile=self.plasma_profile, current_drive=self.current_drive
+            physics=self.physics,
         )
         self.dcll = DCLL(blanket_library=self.blanket_library)
 

process/stellarator.py

@@ -4,7 +4,6 @@
 
 from process.fortran import (
     constants,
-    physics_module as ph,
     stellarator_module as st,
     process_output as po,
     physics_variables,
@@ -25,14 +24,15 @@
     constraint_variables,
     rebco_variables,
     maths_library,
-    profiles_module,
     physics_functions_module,
     neoclassics_module,
     impurity_radiation_module,
+    sctfcoil_module,
 )
 import process.superconductors as superconductors
 import process.physics_functions as physics_funcs
 from process.coolprop_interface import FluidProperties
+from process.physics import rether
 
 logger = logging.getLogger(__name__)
 # Logging handler for console output
@@ -62,6 +62,7 @@ def __init__(
         plasma_profile,
         hcpb,
         current_drive,
+        physics,
     ) -> None:
         """Initialises the Stellarator model's variables
 
@@ -79,6 +80,8 @@ def __init__(
         :type hcpb: process.hcpb.CCFE_HCPB
         :param current_drive: a pointer to the CurrentDrive model, allowing use of CurrentDrives's variables/methods
         :type current_drive: process.current_drive.CurrentDrive
+        :param physics: a pointer to the Physics model, allowing use of Physics's variables/methods
+        :type physics: process.physics.Physics
         """
 
         self.outfile: int = constants.nout
@@ -92,6 +95,7 @@ def __init__(
         self.plasma_profile = plasma_profile
         self.hcpb = hcpb
         self.current_drive = current_drive
+        self.physics = physics
 
     def run(self, output: bool):
         """Routine to call the physics and engineering modules
@@ -110,7 +114,7 @@ def run(self, output: bool):
             self.costs.run()
             self.costs.output()
             self.availability.run(output=True)
-            ph.outplas(self.outfile)
+            self.physics.outplas()
             self.stigma()
             self.stheat(True)
             self.stphys(True)
@@ -203,7 +207,7 @@ def stigma(self):
                 physics_variables.tauei,
                 physics_variables.taueff,
                 physics_variables.powerht,
-            ) = physics_module.pcond(
+            ) = self.physics.pcond(
                 physics_variables.afuel,
                 physics_variables.palpmw,
                 physics_variables.aspect,
@@ -234,7 +238,7 @@ def stigma(self):
                 physics_variables.zeff,
             )
 
-            physics_variables.hfac[iisc] = physics_module.fhfac(i)
+            physics_variables.hfac[iisc] = self.physics.fhfac(i)
 
     def stnewconfig(self):
         """author: J Lion, IPP Greifswald
@@ -1187,12 +1191,6 @@ def stfwbs(self, output: bool):
         volshldo = build_variables.shareaob * build_variables.shldoth
         fwbs_variables.volshld = volshldi + volshldo
 
-        fwbs_variables.whtshld = (
-            fwbs_variables.volshld
-            * fwbs_variables.denstl
-            * (1.0e0 - fwbs_variables.vfshld)
-        )
-
         #  Neutron power lost through holes in first wall (eventually absorbed by
         #  shield)
 
@@ -1728,6 +1726,13 @@ def stfwbs(self, output: bool):
         if fwbs_variables.blktmodel == 0:
             coolvol = coolvol + fwbs_variables.volblkt * fwbs_variables.vfblkt
 
+        # Shield mass
+        fwbs_variables.whtshld = (
+            fwbs_variables.volshld
+            * fwbs_variables.denstl
+            * (1.0e0 - fwbs_variables.vfshld)
+        )
+
         coolvol = coolvol + fwbs_variables.volshld * fwbs_variables.vfshld
 
         #  Penetration shield (set = internal shield)
@@ -2718,6 +2723,7 @@ def stcoil(self, output: bool):
         # NOTE: original implementation used taucq which used a EUROfusion
         # constant in the calculation. This was the minimum allowed quench time.
         # Replacing with the actual quench time.
+        # MN/m^3
         f_vv_actual = (
             2.54e6
             * (3e0 * 1.3e0 * 50e0 * 0.92e0**2e0)
@@ -2736,6 +2742,15 @@ def stcoil(self, output: bool):
             ** (-1)
         )
 
+        # N/m^2
+        # is the vv width the correct length to multiply by to turn the
+        # force density into a stress?
+        sctfcoil_module.vv_stress_quench = (
+            f_vv_actual
+            * 1e6
+            * ((build_variables.d_vv_in + build_variables.d_vv_out) / 2)
+        )
+
         # the conductor fraction is meant of the cable space#
         # This is the old routine which is being replaced for now by the new one below
         #    protect(aio,  tfes,               acs,       aturn,   tdump,  fcond,  fcu,   tba,  tmax   ,ajwpro, vd)
@@ -3250,7 +3265,6 @@ def power_at_ignition_point(self, gyro_frequency_max, te0_available):
         Assumes current peak temperature (which is inaccurate as the cordey pass should be calculated)
         Maybe use this: https://doi.org/10.1088/0029-5515/49/8/085026
         """
-
         te_old = copy(physics_variables.te)
         # Volume averaged physics_variables.te from te0_achievable
         physics_variables.te = te0_available / (1.0e0 + physics_variables.alphat)
@@ -3849,10 +3863,10 @@ def stphys(self, output):
         #  Calculate plasma composition
         # Issue #261 Remove old radiation model
 
-        physics_module.plasma_composition()
+        self.physics.plasma_composition()
 
         # Calculate density and temperature profile quantities
-        profiles_module.plasma_profiles()
+        self.plasma_profile.run()
 
         #  Total field
         physics_variables.btot = np.sqrt(
@@ -4048,7 +4062,7 @@ def stphys(self, output):
 
         #  Calculate ion/electron equilibration power
 
-        physics_variables.piepv = physics_module.rether(
+        physics_variables.piepv = rether(
             physics_variables.alphan,
             physics_variables.alphat,
             physics_variables.dene,
@@ -4174,7 +4188,7 @@ def stphys(self, output):
             physics_variables.tauei,
             physics_variables.taueff,
             physics_variables.powerht,
-        ) = physics_module.pcond(
+        ) = self.physics.pcond(
             physics_variables.afuel,
             physics_variables.palpmw,
             physics_variables.aspect,
@@ -4224,7 +4238,7 @@ def stphys(self, output):
             physics_variables.qfuel,
             physics_variables.rndfuel,
             physics_variables.taup,
-        ) = physics_module.phyaux(
+        ) = self.physics.phyaux(
             physics_variables.aspect,
             physics_variables.dene,
             physics_variables.deni,

tests/unit/test_stellarator.py

@@ -24,6 +24,7 @@
 from process.blanket_library import BlanketLibrary
 from process.fw import Fw
 from process.current_drive import CurrentDrive
+from process.physics import Physics
 
 
 @pytest.fixture
@@ -42,6 +43,7 @@ def stellarator():
         PlasmaProfile(),
         CCFE_HCPB(BlanketLibrary(Fw())),
         CurrentDrive(),
+        Physics(PlasmaProfile(), CurrentDrive()),
     )