Routine that formulates the constraint equations
author: P J Knight (UKAEA)
author: J Morris (UKAEA)
if ieqn is zero or negative, evaluate all the constraint equations, otherwise
evaluate only the ieqnth equation. The code attempts to make for all
equations. All relevant consistency equations should be active in
order to make a self-consistent machine.
References
1.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | m | Number of constraint equations to solve |
||
| integer, | intent(in) | :: | ieqn | Switch for constraint equations to evaluate; |
||
| real(kind=dp), | intent(out), | dimension(m) | :: | cc | Residual error in equation i |
|
| real(kind=dp), | intent(out), | optional | dimension(m) | :: | con | constraint value for equation i in physical units |
| real(kind=dp), | intent(out), | optional | dimension(m) | :: | err | residual error in equation i in physical units |
| character(len=1), | intent(out), | optional | dimension(m) | :: | symbol |
|
| character(len=10), | intent(out), | optional | dimension(m) | :: | units | constraint units in equation i |
subroutine constraint_eqns(m,ieqn,cc,con,err,symbol,units)
!! Routine that formulates the constraint equations
!!
!! **author: P J Knight** (UKAEA)
!!
!! **author: J Morris** (UKAEA)
!!
!! if `ieqn` is zero or negative, evaluate all the constraint equations, otherwise
!! evaluate only the `ieqn`th equation. The code attempts to make \( cc(i) = 0 \) for all
!! \( i \in \{1,\cdots,m\} \) equations. All relevant consistency equations should be active in
!! order to make a self-consistent machine.
!!
!! **References**
!!
!! 1.
use numerics, only: icc
use maths_library, only: variable_error
implicit none
integer, intent(in) :: m
!! Number of constraint equations to solve
integer, intent(in) :: ieqn
!! Switch for constraint equations to evaluate;
real(dp), dimension(m), intent(out) :: cc
!! Residual error in equation i
real(dp), optional, dimension(m), intent(out) :: con
!! constraint value for equation i in physical units
real(dp), optional, dimension(m), intent(out) :: err
!! residual error in equation i in physical units
character(len=1), optional, dimension(m), intent(out) :: symbol
!! `=<`, `>`, `<` symbol for equation i denoting its type
character*10, optional, dimension(m), intent(out) :: units
!! constraint units in equation i
! Local variables
integer :: i, i1, i2
real(dp) :: tmp_cc = 0
!! Residual error in constraint equation
real(dp) :: tmp_con = 0
!! constraint value for constraint equation in physical units
real(dp) :: tmp_err = 0
!! residual error in constraint equation in physical units
character(len=1) :: tmp_symbol = ' '
!! `=<`, `>`, `<` symbol for constraint equation denoting its type
character(len=10) :: tmp_units = ' '
!! constraint units in constraint equation
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! If ieqn is positive, only evaluate the 'ieqn'th constraint residue,
! otherwise evaluate all m constraint residues
if (ieqn > 0) then
i1 = ieqn ; i2 = ieqn
else
i1 = 1 ; i2 = m
end if
! Consistency (equality) constraints should converge to zero.
do i = i1,i2
! The constraint value in physical units is
! a) for consistency equations, the quantity to be equated, or
! b) for limit equations, the limiting value.
! The symbol is = for a consistency equation, < for an upper limit
! or > for a lower limit.
select case (icc(i))
! Relationship between beta, temperature (keV) and density
case (1); call constraint_eqn_001(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Global plasma power balance equation
case (2); call constraint_eqn_002(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Global power balance equation for ions
case (3); call constraint_eqn_003(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Global power balance equation for electrons
case (4); call constraint_eqn_004(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for density upper limit
case (5); call constraint_eqn_005(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for epsilon beta-poloidal upper limit
case (6); call constraint_eqn_006(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for hot beam ion density
case (7); call constraint_eqn_007(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for neutron wall load upper limit
case (8); call constraint_eqn_008(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for fusion power upper limit
case (9); call constraint_eqn_009(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Obsolete
case (10); call constraint_eqn_010(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for radial build
case (11); call constraint_eqn_011(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for volt-second capability lower limit
case (12); call constraint_eqn_012(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for burn time lower limit
case (13); call constraint_eqn_013(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation to fix number of NBI decay lengths to plasma centre
case (14); call constraint_eqn_014(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for L-H power threshold limit
case (15); call constraint_eqn_015(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for net electric power lower limit
case (16); call constraint_eqn_016(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for radiation power upper limit
case (17); call constraint_eqn_017(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for divertor heat load upper limit
case (18); call constraint_eqn_018(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for MVA upper limit
case (19); call constraint_eqn_019(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for neutral beam tangency radius upper limit
case (20); call constraint_eqn_020(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for minor radius lower limit
case (21); call constraint_eqn_021(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for divertor collision/connection length ratio upper limit
case (22); call constraint_eqn_022(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for conducting shell radius / rminor upper limit
case (23); call constraint_eqn_023(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for beta upper limit
case (24); call constraint_eqn_024(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for peak toroidal field upper limit
case (25); call constraint_eqn_025(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for Central Solenoid current density upper limit at EOF
case (26); call constraint_eqn_026(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for Central Solenoid current density upper limit at BOP
case (27); call constraint_eqn_027(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for fusion gain (big Q) lower limit
case (28); call constraint_eqn_028(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for inboard major radius
case (29); call constraint_eqn_029(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for injection power upper limit
case (30); call constraint_eqn_030(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for TF coil case stress upper limit (SCTF)
case (31); call constraint_eqn_031(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for TF coil conduit stress upper limit (SCTF)
case (32); call constraint_eqn_032(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for TF coil operating/critical J upper limit (SCTF)
case (33); call constraint_eqn_033(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for TF coil dump voltage upper limit (SCTF)
case (34); call constraint_eqn_034(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for TF coil J_wp/J_prot upper limit (SCTF)
case (35); call constraint_eqn_035(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for TF coil s/c temperature margin lower limit (SCTF)
case (36); call constraint_eqn_036(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for current drive gamma upper limit
case (37); call constraint_eqn_037(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for first wall temperature upper limit
case (39); call constraint_eqn_039(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for auxiliary power lower limit
case (40); call constraint_eqn_040(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for plasma current ramp-up time lower limit
case (41); call constraint_eqn_041(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for cycle time lower limit
case (42); call constraint_eqn_042(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for average centrepost temperature
case (43); call constraint_eqn_043(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for centrepost temperature upper limit (TART)
case (44); call constraint_eqn_044(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for edge safety factor lower limit (TART)
case (45); call constraint_eqn_045(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for Ip/Irod upper limit (TART)
case (46); call constraint_eqn_046(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for TF coil toroidal thickness upper limit
case (47); call constraint_eqn_047(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for poloidal beta upper limit
case (48); call constraint_eqn_048(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Issue #508 Remove RFP option Equation to ensure reversal parameter F is negative
case (49); call constraint_eqn_049(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! IFE option: Equation for repetition rate upper limit
case (50); call constraint_eqn_050(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation to enforce startup flux = available startup flux
case (51); call constraint_eqn_051(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for tritium breeding ratio lower limit
case (52); call constraint_eqn_052(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for fast neutron fluence on TF coil upper limit
case (53); call constraint_eqn_053(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for peak TF coil nuclear heating upper limit
case (54); call constraint_eqn_054(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for helium concentration in vacuum vessel upper limit
case (55); call constraint_eqn_055(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for power through separatrix / major radius upper limit
case (56); call constraint_eqn_056(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Obsolete
case (57); call constraint_eqn_057(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Obsolete
case (58); call constraint_eqn_058(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for neutral beam shine-through fraction upper limit
case (59); call constraint_eqn_059(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for Central Solenoid s/c temperature margin lower limit
case (60); call constraint_eqn_060(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for availability limit
case (61); call constraint_eqn_061(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Lower limit on taup/taueff the ratio of alpha particle to energy confinement times
case (62); call constraint_eqn_062(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Upper limit on niterpump (vacuum_model = simple)
case (63); call constraint_eqn_063(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Upper limit on Zeff
case (64); call constraint_eqn_064(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Limit TF dump time to calculated quench time
case (65); call constraint_eqn_065(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Limit on rate of change of energy in poloidal field
case (66); call constraint_eqn_066(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Simple upper limit on radiation wall load
case (67); call constraint_eqn_067(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! New Psep scaling (PsepB/qAR)
case (68); call constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Ensure separatrix power is less than value from Kallenbach divertor
case (69); call constraint_eqn_069(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Separatrix temperature consistency
case (70); call constraint_eqn_070(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Separatrix density consistency
case (71); call constraint_eqn_071(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Central Solenoid Tresca yield criterion
case (72); call constraint_eqn_072(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! ensure separatrix power is greater than the L-H power + auxiliary power
case (73); call constraint_eqn_073(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! ensure TF coil quench temperature < tmax_croco
case (74); call constraint_eqn_074(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! ensure that TF coil current / copper area < Maximum value ONLY used for croco HTS coil
case (75); call constraint_eqn_075(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Eich critical separatrix density model
case (76); call constraint_eqn_076(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for maximum TF current per turn upper limit
case (77); call constraint_eqn_077(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for Reinke criterion, divertor impurity fraction lower limit
case (78); call constraint_eqn_078(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Equation for maximum CS field
case (79); call constraint_eqn_079(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Lower limit pdivt
case (80); call constraint_eqn_080(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint equation making sure that ne(0) > ne(ped)
case (81); call constraint_eqn_081(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint equation making sure that stellarator coils dont touch in toroidal direction
case (82); call constraint_eqn_082(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint ensuring radial build consistency for stellarators
case (83); call constraint_eqn_083(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint for lower limit of beta
case (84); call constraint_eqn_084(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint for CP lifetime
case (85); call constraint_eqn_085(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint for turn dimension
case (86); call constraint_eqn_086(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint for cryogenic power
case (87); call constraint_eqn_087(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint for TF coil strain
case (88); call constraint_eqn_088(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! ensure that OH coil current / copper area < Maximum value ONLY used for croco HTS coil
case (89); call constraint_eqn_089(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint for minimum CS stress load cycles
case (90); call constraint_eqn_090(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
! Constraint for indication of ECRH ignitability
case (91); call constraint_eqn_091(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
case default
idiags(1) = icc(i)
call report_error(13)
tmp_cc = 0
tmp_con = 0
tmp_err = 0
tmp_symbol = ' '
tmp_units = ' '
end select
cc(i) = tmp_cc
if (present(con)) then
con(i) = tmp_con
if (present(err)) err(i) = tmp_err
if (present(symbol)) symbol(i) = tmp_symbol
if (present(units)) units(i) = tmp_units
end if
! Crude method of catching NaN errors
!if ((abs(cc(i)) > 9.99D99).or.(cc(i) /= cc(i))) then
if (variable_error(cc(i))) then
! Add debugging lines as appropriate...
select case (icc(i))
! Relationship between beta, temperature (keV) and density (consistency equation)
case (1); call constraint_err_001()
! Equation for net electric power lower limit
case (16); call constraint_err_016()
! Equation for injection power upper limit
case (30); call constraint_err_030()
! Limit on rate of change of energy in poloidal field
case (66); call constraint_err_066()
end select
idiags(1) = icc(i) ; fdiags(1) = cc(i)
call report_error(14)
end if
end do
! Issue 505 Reverse the sign so it works as an inequality constraint (cc(i) > 0)
! This will have no effect if it is used as an equality constraint because it will be squared.
cc = -cc
end subroutine constraint_eqns