Initial commit from Polytechnique Montreal

1 files changed, 356 insertions, 0 deletions

diff --git a/Donjon/data/rep900_sim_recopy_proc/Steady.c2m b/Donjon/data/rep900_sim_recopy_proc/Steady.c2m
new file mode 100644
index 0000000..f21903d
--- /dev/null
+++ b/Donjon/data/rep900_sim_recopy_proc/Steady.c2m
@@ -0,0 +1,356 @@
+*****************************************************************
+*                                                               *
+* Procedure :  Steady.c2m                                       *
+* Purpose   :  reactor steady state computation including       *
+*              neutron physics, thermalhydraulics, local xenon  *
+*              effects and critical boron search feature        *
+* Author    :  V. Salino                                        *
+*                                                               *
+* CALL      :                                                   *
+*  Power Fmap Matex Thermo MicroF Burn := Steady                *
+*        Fmap Matex Thermo MicroF Burn ReflXS Track ::          *
+*        <<CritCB>> <<CBinit>>                                  *
+*        >>CB<< >>keff<< ;                                      *
+*                                                               *
+* Parameters: input                                             *
+*  CritCB    set to True to search critical boron concentration *
+*            computation, False otherwise (logical variable).   *
+*  CBinit    initialization for boron concentration (ppm) when  *
+*            CritCB.EQ.True, otherwise boron concentration to   *
+*            be used.                                           *
+*  Prel      relative power=requested power/nominal power.      *
+*  InletTmod inlet coolant temperature (Celsius).               *
+*                                                               *
+* Parameters: output                                            *
+*  CB        critical boron concentration if CritCB.EQ.True,    *
+*            can be negative in which case CB is extrapolated   *
+*            and should not considered physical but fictive.    *
+*  keff      k-effective, differs from 1 only if CB is          *
+*            negative or if CritCB.EQ.False.                    *
+*                                                               *
+*****************************************************************
+PARAMETER Power Flux Thermo MicroF Burn Fmap Matex SapUOX SapMOX 
+          MacroRefl Track ::
+  ::: LINKED_LIST Power Flux Thermo MicroF Burn Fmap Matex SapUOX 
+                  SapMOX MacroRefl Track ; ;
+MODULE SCR: NCR: MACINI: RESINI: TRIVAA: FLUD: FLPOW: THM: EVO: SIM:
+       GREP: DELETE: ABORT: END: ;
+*--
+* Local objects and procedure
+*--
+LINKED_LIST FluxAv MacroTot MacroF System PowerIter 
+            OldPower
+             ;
+PROCEDURE ComparePdist ;
+*--
+* Procedure arguments
+*--
+LOGICAL CritCB evo pred ;
+INTEGER BUindex ;
+REAL CB ;
+  :: >>CritCB<< >>CB<< >>evo<< >>pred<< ;
+*--
+* Local variables
+*--
+REAL CBmax := 2000.0 ;
+REAL CBmin :=    0.0 ;
+REAL ErrorRho ErrorCB ErrorTFuel ErrorDCool ErrorTCool ErrorPdis ;
+REAL EpsRho EpsCB EpsTFuel EpsDCool EpsPdis :=
+     0.5    0.05  0.5      5.0E-5   5.0E-4   ;
+INTEGER Iter := 1 ;
+REAL CBinterp CBp1 CBm1 keff DeltaRho Rho Rhom1 ;
+LOGICAL CONV := $False_L ;
+REAL DiffBorWorth ;
+*--
+* Thermalhydraulics parameters
+*--
+REAL dx       := 21.5 ;
+REAL Tot_tub  := 6.6E-03 2.0 ** $Pi_R * 25.0  * ;
+REAL Tot_pin  := 4.7E-03 2.0 ** $Pi_R * 264.0 * ;
+REAL asssect  := dx dx * 1.E-04 * Tot_tub - Tot_pin - ;
+REAL coresect := 157.0 asssect * ;
+REAL sass     := dx dx * 1.E-04 * ;
+*--
+* Initialized with a THM computation using a flat power distribution
+*--
+REAL Ptot := 2750.0 ;
+REAL BundPow := Ptot 157.0 / 29.0 / 1.0E+3 * (* MW to kW *) ;
+Fmap := RESINI: Fmap ::
+  BUNDLE-POW SAME <<BundPow>> ;
+*--
+* Iteration loop with feedback, flux, xenon and simplified
+* thermalhydraulics
+*--
+IF pred THEN
+FluxAv := Flux ;
+ENDIF ;
+*--
+REPEAT
+*--
+* Determine CB to use for interpolation
+*--
+  IF CB CBmin < THEN
+    EVALUATE CBinterp := CBmin ;
+  ELSEIF CB CBmax > THEN
+    EVALUATE CBinterp := CBmax ;
+  ELSE
+    EVALUATE CBinterp := CB ;
+  ENDIF ;
+*--
+* Thermalhydraulics computation, xenon saturation
+*--
+  IF Iter 1 = THEN
+    Fmap := SIM: Fmap ::
+    EDIT 1
+    SET-PARAM 'C-BORE' <<CBinterp>> ;
+  ELSE
+  System := DELETE: System ;
+    Thermo Fmap := THM: Thermo Fmap :: EDIT 0
+      SET-PARAM 'C-BORE' <<CBinterp>> ;
+    ! Evolution jusqu'a l'equilibre seulement si demande
+    ! Sinon on va rechercher la CB crit sans etre a l'equilibre Xe
+    IF evo THEN
+      ECHO "bugg1" ;
+      Burn MicroF := EVO: Burn MicroF PowerIter ::
+        EDIT 0 FLUX_POW DEPL 100. DAY KEEP ;
+      ECHO "bugg2" ;
+      Burn MicroF := EVO: Burn MicroF PowerIter ::
+        EDIT 0 FLUX_POW SAVE 100. DAY KEEP ;
+    ENDIF ;
+  ENDIF ;
+*--
+* Saphyb Interpolation
+* Historic effects are neglected for simplicity's sake
+*--
+  MicroF := SCR: MicroF SapUOX SapMOX Fmap ::
+      EDIT 0 RES
+      MICRO LINEAR
+      TABLE SapUOX 'BURN'
+      MIX 4
+         INST-BURN
+        SET TF  526.85 ! Celsius
+        SET TCA 326.85 ! Celsius
+        SET DCA 0.659 ! g/cm3
+        SET ppmB 500.0 ! ppm
+        ADD 'TF' 526.85 MAP
+          REF 'BURN' SAMEASREF ENDREF
+        ADD 'TCA' 326.85 MAP
+          REF 'BURN' SAMEASREF ENDREF
+        ADD 'DCA' 0.659 MAP
+          REF 'BURN' SAMEASREF ENDREF
+        ADD 'ppmB' 500.0 MAP
+          REF 'BURN' SAMEASREF ENDREF
+        MICRO ALL
+* Recherche poids Xenon : commenter/decommenter
+!             'Xe135'   *
+!             'I135'    0.
+              'B10'     *
+              'B11'     *
+* Temporairement, Sm/Pm sont remis a leurs valeurs de Saphyb
+* (equivalent a ne pas les particulariser). Apres premier
+* appel a Steady, on pourrait enlever cela pour les laisser
+* evoluer librement.
+!             'Nd147'   0.
+!             'Pm147'   0.
+!             'Pm148'   0.
+!             'Pm148m'  0.
+!             'Pm149'   0.
+!             'Sm149'   0.
+             '*MAC*RES' 1.
+      ENDMIX
+      TABLE SapMOX 'BURN'
+      MIX 5
+         INST-BURN
+        SET TF  526.85 ! Celsius
+        SET TCA 326.85 ! Celsius
+        SET DCA 0.659 ! g/cm3
+        SET ppmB 500.0 ! ppm
+        ADD 'TF' 526.85 MAP
+          REF 'BURN' SAMEASREF ENDREF
+        ADD 'TCA' 326.85 MAP
+          REF 'BURN' SAMEASREF ENDREF
+        ADD 'DCA' 0.659 MAP
+          REF 'BURN' SAMEASREF ENDREF
+        ADD 'ppmB' 500.0 MAP
+          REF 'BURN' SAMEASREF ENDREF
+        MICRO ALL
+* Recherche poids Xenon : commenter/decommenter
+!             'Xe135'   *
+!             'I135'    0.
+              'B10'     *
+              'B11'     *
+* Temporairement, Sm/Pm sont remis a leurs valeurs de Saphyb
+* (equivalent a ne pas les particulariser). Apres premier
+* appel a Steady, on pourrait enlever cela pour les laisser
+* evoluer librement.
+!             'Nd147'   0.
+!             'Pm147'   0.
+!             'Pm148'   0.
+!             'Pm148m'  0.
+!             'Pm149'   0.
+!             'Sm149'   0.
+             '*MAC*RES' 1.
+      ENDMIX
+      CHAIN
+        B10      STABLE
+        B11      STABLE
+        I135     NG    0.
+        Xe135    NG    0. FROM DECAY 1.0E+00 I135
+        Nd147   STABLE
+        Pm147   STABLE
+        Pm148   STABLE
+        Pm148m  STABLE
+        Pm149   STABLE
+        Sm149   STABLE
+*--------- With Samarium -----------------------------------
+             Nd147    NG    0.
+             Pm147    NG    0. FROM DECAY 1.0E+00 Nd147
+             Pm148    NG    0. FROM NG    5.3E-01 Pm147
+             Pm148m   NG    0. FROM NG    4.7E-01 Pm147
+             Pm149    NG    0. FROM NG    1.0E+00 Pm148
+                                    NG    1.0E+00 Pm148m
+             Sm149    NG    0. FROM DECAY 1.0E+00 Pm149
+        MACT     NFTOT 0.
+      ENDCHAIN
+  ;
+  MacroF := MicroF :: STEP UP "MACROLIB" ;
+*--
+* Flux and power computation
+*--
+  MacroTot Matex := MACINI: Matex MacroRefl MacroF ;
+  MacroF := DELETE: MacroF ;
+  System := TRIVAA: MacroTot Track :: EDIT 0 ;
+  MacroTot := DELETE: MacroTot ;
+  IF Iter 1 = THEN
+    Flux := FLUD: Flux System Track :: EDIT 0
+      ADI 4 EXTE 1.0E-5 ACCE 5 3 ;
+  ELSE
+    Flux := FLUD: Flux System Track :: EDIT 0
+      RELAX 0.15 ;
+  ENDIF ;
+  GREP: Flux :: GETVAL 'K-EFFECTIVE ' 1  >>keff<<  ;
+  EVALUATE Rho := 1. 1. keff / - 1.0E5 * ;
+  IF Iter 1 > THEN
+    PowerIter := DELETE: PowerIter ;
+  ENDIF ;
+  PowerIter Fmap := FLPOW: Fmap Flux Track Matex ::
+    EDIT 0 PTOT <<Ptot>> PRINT DISTR POWER ;
+*--
+* Compute convergence errors (absolute values) on :
+* - reactivity (pcm)
+* - boron concentration (ppm)
+* - 3D power distribution (relative discrepancy, maximum)
+* - fuel temperature distribution (K, maximum)
+* - moderator density distribution (g/cm3, maximum)
+*--
+  IF Iter 1 > THEN
+    EVALUATE ErrorRho := Rho Rhom1 - ABS ;
+    EVALUATE ErrorCB  := CB CBm1 - ABS ;
+    Fmap OldPower PowerIter := ComparePdist :: >>ErrorPdis<< ;
+    OldPower := DELETE: OldPower ;
+    GREP: Thermo :: GETVAL 'ERROR-T-FUEL' 1  >>ErrorTFuel<<  ;
+    GREP: Thermo :: GETVAL 'ERROR-D-COOL' 1  >>ErrorDCool<<  ;
+  ENDIF ;
+  OldPower := PowerIter ;
+*--
+* Reestimate critical CB. 3rd iteration is waited, in order to have a
+* not-too-false k-effective, in particular regarding thermalhydraulics
+* and neutron physics convergence. Otherwise predicted CB would be too
+* far from the truth.
+*--
+  IF CritCB Iter 2 > * THEN
+* To evaluate DiffBorWorth, we have 3 other possibilities :
+* - Evaluate properly DiffBorWorth (internship ?), may be costly
+* - Use the 2 previous computations to find (1 line with 2 points)
+*   (may be instable during the first iterations but effective
+*   near the end)
+* - Use module FIND0 (may be long)
+* Plot reactivity=f(CB), see which is going to be good
+    EVALUATE DiffBorWorth := -7.0 ; ! pcm/ppm
+    EVALUATE CBp1 := CBinterp Rho DiffBorWorth / - ;
+  ENDIF ;
+*--
+* Prints
+*--
+  ECHO "#---" ;
+  ECHO "Iter" Iter ;
+  ECHO "keff=" keff "CB=" CB "Rho=" Rho ;
+  ECHO "CBinterp=" CBinterp ;
+  IF CritCB Iter 2 > * THEN
+    ECHO "DiffBorWorth=" DiffBorWorth "CBp1" CBp1 ;
+  ENDIF ;
+  ECHO "---" ;
+*--
+* Check convergence
+*--
+  IF Iter 3 > THEN
+    ECHO "Discrepancies between the previous iteration and current"
+         "iteration" ;
+    ECHO " Rho (pcm)     CB (ppm)      TFuel (K)     DCool (K)    "
+         "Pdis (relat)" ;
+    ECHO ErrorRho ErrorCB ErrorTFuel ErrorDCool ErrorPdis ;
+    ECHO "" ;
+    ECHO "Convergence criterion (epsilon)" ;
+    ECHO " Rho (pcm)     CB (ppm)      TFuel (K)     DCool (K)    "
+         "Pdis (relat)" ;
+    ECHO EpsRho EpsCB EpsTFuel EpsDCool EpsPdis ;
+    ECHO "" ;
+    ECHO "Convergence reached ?" ;
+    ECHO " Rho (pcm)     CB (ppm)      TFuel (K)     DCool (K)    "
+         "Pdis (relat)" ;
+    ECHO "" ErrorRho EpsRho < "      " ErrorCB EpsCB < "      "
+         ErrorTFuel EpsTFuel < "       " ErrorDCool EpsDCool <
+         "      " ErrorPdis  EpsPdis  < ;
+    EVALUATE CONV := ErrorRho   EpsRho   <
+                     ErrorCB    EpsCB    < *
+                     ErrorTFuel EpsTFuel < *
+                     ErrorDCool EpsDCool < *
+                     ErrorPdis  EpsPdis  < * ;
+  ENDIF ;
+*--
+* Prepare next iteration
+*--
+  IF CONV NOT THEN
+    EVALUATE CBm1 Rhom1 := CB Rho ;
+    IF CritCB Iter 2 > * THEN
+      EVALUATE CB := CBp1 ;
+    ENDIF ;
+    EVALUATE Iter := Iter 1 + ;
+  ENDIF ;
+*--
+  IF Iter 100 > THEN
+    ECHO "Steady.c2m: maximum iteration reached (50)." ;
+    ABORT: ;
+  ENDIF ;
+UNTIL Iter 4 = ;
+*--
+* Predicteur-correcteur
+*--
+IF pred THEN
+  PowerIter := DELETE: PowerIter ;
+  FluxAv := FLUD: FluxAv System Track :: EDIT 0
+      RELAX 0.5 ;
+  PowerIter Fmap := FLPOW: Fmap FluxAv Track Matex ::
+    EDIT 0 PTOT <<Ptot>> PRINT DISTR POWER ;
+  FluxAv := DELETE: FluxAv ;
+ENDIF ;
+
+*--
+* Print Xe and I average fission yields. This average depends
+* on burnup, thermalhydraulics and xenon conditions.
+*--
+REAL YieldI YieldXe ;
+GREP: MicroF :: STEP UP 'DEPL-CHAIN'
+                   GETVAL 'FISSIONYIELD' 1 >>YieldI<<
+                   GETVAL 'FISSIONYIELD' 2 >>YieldXe<< ;
+ECHO "Mean fission yield for I =" YieldI  ;
+ECHO "Mean fission yield for Xe=" YieldXe ;
+
+Power := PowerIter ;
+PowerIter := DELETE: PowerIter ;
+System := DELETE: System ;
+
+  :: <<CB>> ;
+  :: <<keff>> ; 
+
+END: ;