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 cc(i)=0 for all
i∈{1,⋯,m} 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
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)
! Obsolete
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 f_alpha_energy_confinement 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)
! Constraint for D/T ratio
case (92); call constraint_eqn_092(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
if (isnan(cc(i)).or.(abs(cc(i))>9.99D99)) 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