Difference between revisions of "MSCS for burnup calculations"

From Serpent Wiki
Jump to: navigation, search
(→‎Description)
(→‎Setup)
 
(19 intermediate revisions by 3 users not shown)
Line 13: Line 13:
 
One of the provided inputs uses the explicit Euler discretization for the depletion solution with stochastic approximation based relaxation applied to the tallied transmutation cross sections and one-group fluxes. The second input uses the Stochastic Implicit Euler's method<ref>J. Dufek and H. Anglart, "''Derivation of a stable coupling scheme for Monte Carlo burnup calculations with the thermal-hydraulic feedback''", Ann. Nucl. Energy '''62''', pp. 260 (2013) DOI: [http://dx.doi.org/10.1016/j.anucene.2013.06.025 10.1016/j.anucene.2013.06.025]</ref>.
 
One of the provided inputs uses the explicit Euler discretization for the depletion solution with stochastic approximation based relaxation applied to the tallied transmutation cross sections and one-group fluxes. The second input uses the Stochastic Implicit Euler's method<ref>J. Dufek and H. Anglart, "''Derivation of a stable coupling scheme for Monte Carlo burnup calculations with the thermal-hydraulic feedback''", Ann. Nucl. Energy '''62''', pp. 260 (2013) DOI: [http://dx.doi.org/10.1016/j.anucene.2013.06.025 10.1016/j.anucene.2013.06.025]</ref>.
  
−
== Files ==
+
== Files for explicit Euler ==
  
−
=== MSCS.py ===
+
=== MSCS.py (explicit Euler) ===
  
 
  <nowiki>#############################################################
 
  <nowiki>#############################################################
 
#                                                          #
 
#                                                          #
−
#          Minimal Serpent Coupling Script v 0.2           #
+
#          Minimal Serpent Coupling Script v 0.3           #
 
#          For burnup example                              #
 
#          For burnup example                              #
 
#                                                          #
 
#                                                          #
 
# Created by:  Ville Valtavirta              2016/01/04    #
 
# Created by:  Ville Valtavirta              2016/01/04    #
−
# Modified by: Ville Valtavirta              2016/01/05     #
+
# Modified by: Ville Valtavirta              2018/10/26     #
 
#                                                          #
 
#                                                          #
 
#############################################################
 
#############################################################
Line 34: Line 34:
 
# Path to Serpent executable
 
# Path to Serpent executable
  
−
sssexe = '/home/vvvillehe/Serpent2/2.1.27/sss2'
+
sssexe = '/home/vvvillehe/Serpent2/2.1.31/sss2'
  
 
#######################################################
 
#######################################################
Line 211: Line 211:
  
 
             if int(line) != -1:
 
             if int(line) != -1:
−
                 if int(line) == signal.SIGUSR1:
+
                 if int(line) == signal.SIGUSR1.value:
 
                     # Got the signal to resume
 
                     # Got the signal to resume
  
 
                     sleeping = 0
 
                     sleeping = 0
  
−
                 elif int(line) == signal.SIGUSR2:
+
                 elif int(line) == signal.SIGUSR2.value:
 
                     # Got the signal to move to next time point
 
                     # Got the signal to move to next time point
  
 
                     iterating = 0
 
                     iterating = 0
 
                     sleeping = 0
 
                     sleeping = 0
−
                 elif int(line) == signal.SIGTERM:
+
                 elif int(line) == signal.SIGTERM.value:
 
                     # Got the signal to end the calculation
 
                     # Got the signal to end the calculation
  
Line 230: Line 230:
 
                     # Unknown signal
 
                     # Unknown signal
  
−
                     print "\nUnknown singal read from file, exiting\n"
+
                     print("\nUnknown signal read from file, exiting\n")
  
 
                     # Exit
 
                     # Exit
Line 266: Line 266:
  
 
             P.append(float(strtuple[8]))
 
             P.append(float(strtuple[8]))
 +
 +
        file_in.close()
  
 
         if iterating == 0:
 
         if iterating == 0:
Line 330: Line 332:
 
         #########################################
 
         #########################################
  
−
         print "\nMSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS"
+
         print("\nMSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS")
−
         print "|                                          |"
+
         print("|                                          |")
−
         print "|  Solving fuel temperature for t = {:1d} d  |".format(curdays)
+
         print("|  Solving fuel temperature for t = {:1d} d  |".format(curdays))
−
         print "|                                          |"
+
         print("|                                          |")
−
         print "MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS"
+
         print("MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS")
  
 
         for i in range(10):
 
         for i in range(10):
Line 397: Line 399:
 
         file_out = open('./com.in','w')
 
         file_out = open('./com.in','w')
  
−
         file_out.write(str(signal.SIGUSR1))
+
         file_out.write(str(signal.SIGUSR1.value))
  
 
         file_out.close()
 
         file_out.close()
Line 414: Line 416:
 
     file_out = open('./com.in','w')
 
     file_out = open('./com.in','w')
  
−
     file_out.write(str(signal.SIGUSR1))
+
     file_out.write(str(signal.SIGUSR1.value))
  
 
     file_out.close()</nowiki>
 
     file_out.close()</nowiki>
Line 431: Line 433:
 
%    proceeds otherwise normally. This can be used to check that
 
%    proceeds otherwise normally. This can be used to check that
 
%    the calculation sequence itself runs to completion
 
%    the calculation sequence itself runs to completion
−
set acelib "/home/vvvillehe/XSdata/xsdir0K"
+
 
 
%set declib "<path-here>"
 
%set declib "<path-here>"
 
%set nfylib "<path-here>"
 
%set nfylib "<path-here>"
Line 466: Line 468:
 
% --- Axial depletion zone division for the fuel material
 
% --- Axial depletion zone division for the fuel material
  
−
div fuel subz 10 -200.0 200.0
+
div fuel subz 10 -100.0 100.0
  
 
% --- Cladding material:
 
% --- Cladding material:
Line 551: Line 553:
 
fuel  8  11.808
 
fuel  8  11.808
 
fuel  9  11.808
 
fuel  9  11.808
−
fuel 10  11.808</nowiki>
+
fuel 10  11.808
 +
 
 +
% --- Disable group constant generation
 +
 
 +
set gcu -1</nowiki>
 +
 
 +
== Files for stochastic implicit Euler ==
 +
 
 +
=== MSCS.py (Stochastic Implicit Euler) ===
 +
 
 +
<nowiki>#############################################################
 +
#                                                          #
 +
#          Minimal Serpent Coupling Script v 0.2            #
 +
#          For burnup example (Stochastic Implicit Euler)  #
 +
#                                                          #
 +
# Created by:  Ville Valtavirta              2016/01/04    #
 +
# Modified by: Ville Valtavirta              2020/10/26    #
 +
#                                                          #
 +
#############################################################
 +
 
 +
import os
 +
import signal
 +
import math
 +
import time
 +
 
 +
# Path to Serpent executable
 +
 
 +
sssexe = '/home/vvvillehe/Serpent2/2.1.31/sss2'
 +
 
 +
#######################################################
 +
# Create the Serpent input-file for this run          #
 +
# (process id or communication file must be appended) #
 +
#######################################################
 +
 
 +
# Open original input for reading
 +
 
 +
file_in = open('./input','r')
 +
 
 +
# Open a new input file for writing
 +
 
 +
file_out = open('./coupledinput','w')
 +
 
 +
# Write original input to new file
 +
 
 +
for line in file_in:
 +
    file_out.write(line)
 +
 
 +
# Close original input file
 +
 
 +
file_in.close()
 +
 
 +
# Append signalling mode
 +
 
 +
file_out.write('\n')
 +
file_out.write('set comfile com.in com.out\n')
 +
 
 +
# Append interface names
 +
 
 +
file_out.write('\n')
 +
file_out.write('ifc cool.ifc\n\n')
 +
file_out.write('ifc fuel.ifc\n')
 +
 
 +
# Close new input file
 +
 
 +
file_out.close()
 +
 
 +
##############################################
 +
# Write the initial coolant interface file  #
 +
# (Coolant conditions will be held constant) #
 +
##############################################
 +
 
 +
file_out = open('./cool.ifc','w')
 +
 
 +
# Write the header line (TYPE MAT OUT)
 +
 
 +
file_out.write('2 cool 1\n')
 +
 
 +
# Write the output line (OUTFILE NZ ZMIN ZMAX NR)
 +
 
 +
file_out.write('coolifc.out 10 -100 100 1\n')
 +
 
 +
# Write the mesh type
 +
 
 +
file_out.write('1\n')
 +
 
 +
# Write the mesh size (NX XMIN XMAX NY YMIN YMAX NZ ZMIN ZMAX)
 +
 
 +
file_out.write('1 -0.75 0.75 1 -0.75 0.75 1 -100 100\n')
 +
 
 +
# Write initial coolant temperatures and densities
 +
 
 +
for i in range(1):
 +
    file_out.write('-0.9  520.0\n')
 +
 
 +
# Close interface file
 +
 
 +
file_out.close()
 +
 
 +
##############################################
 +
# Write the initial fuel interface file      #
 +
# (Fuel temperature will be updated)        #
 +
##############################################
 +
 
 +
file_out = open('./fuel.ifc','w')
 +
 
 +
# Write the header line (TYPE MAT OUT)
 +
 
 +
file_out.write('2 fuel 0\n')
 +
 
 +
# Write the mesh type
 +
 
 +
file_out.write('1\n')
 +
 
 +
# Write the mesh size (NX XMIN XMAX NY YMIN YMAX NZ ZMIN ZMAX)
 +
 
 +
file_out.write('1 -0.75 0.75 1 -0.75 0.75 10 -100 100\n')
 +
 
 +
# Write initial fuel temperatures and densities
 +
 
 +
for i in range(10):
 +
    file_out.write('-10.424 310.0\n')
 +
 
 +
# Close interface file
 +
 
 +
file_out.close()
 +
 
 +
# Archive the initial fuel interface
 +
 
 +
os.system('cp ./fuel.ifc ./fuel.ifc0')
 +
 
 +
################################
 +
# Start the Serpent simulation #
 +
################################
 +
 
 +
# Create a command string that will start the Serpent simulation
 +
 
 +
runcommand = sssexe+' -omp 3 ./coupledinput &'
 +
 
 +
# Execute the command string
 +
 
 +
os.system(runcommand)
 +
 
 +
#########################################
 +
# Initialize the fuel behavior solution #
 +
#########################################
 +
 
 +
TFU = []
 +
BU  = []
 +
 
 +
for i in range(10):
 +
    TFU.append(300.0)
 +
    BU.append(0.0)
 +
 
 +
# Reset time step
 +
 
 +
curtime = 0
 +
curdays = 0
 +
 
 +
# Depletion step lengths (in days)
 +
 
 +
steplengths = [2, 3]
 +
 
 +
#############################
 +
# Loop over depletion steps #
 +
#############################
 +
 
 +
simulating = 1
 +
 
 +
while simulating == 1:
 +
 
 +
    ####################################################
 +
    # Picard iteration loop for current depletion step #
 +
    ####################################################
 +
 
 +
    iterating = 1
 +
 
 +
    while iterating == 1:
 +
        ###################
 +
        # Wait for signal #
 +
        ###################
 +
 
 +
        sleeping = 1
 +
 
 +
        while sleeping == 1:
 +
 
 +
            # Sleep for two seconds
 +
 
 +
            time.sleep(2)
 +
 
 +
            # Open file to check if we got a signal
 +
 
 +
            fin = open('./com.out','r')
 +
 
 +
            # Read line
 +
 
 +
            line = fin.readline()
 +
 
 +
            # Close file
 +
 
 +
            fin.close()
 +
 
 +
            # Check signal
 +
 
 +
            if int(line) != -1:
 +
                if int(line) == signal.SIGUSR1.value:
 +
                    # Got the signal to resume
 +
 
 +
                    sleeping = 0
 +
 
 +
                elif int(line) == signal.SIGUSR2.value:
 +
                    # Got the signal to move to next time point
 +
 
 +
                    iterating = 0
 +
                    sleeping = 0
 +
                elif int(line) == signal.SIGTERM.value:
 +
                    # Got the signal to end the calculation
 +
 
 +
                    iterating = 0
 +
                    sleeping = 0
 +
                    simulating = 0
 +
                else:
 +
                    # Unknown signal
 +
 
 +
                    print("\nUnknown signal read from file, exiting\n")
 +
 
 +
                    # Exit
 +
 
 +
                    quit()
 +
 
 +
                # Reset the signal in the file
 +
 
 +
                file_out = open('./com.out','w')
 +
 
 +
                file_out.write('-1')
 +
 
 +
                file_out.close()
 +
 
 +
        ###########################
 +
        # Read power distribution #
 +
        ###########################
 +
 
 +
        # Reset power distribution
 +
 
 +
        P = []
 +
 
 +
        # Get output file number corresponding for current time
 +
        # Serpent versions up to 2.1.29 use the following numbering with SIE:
 +
        # Predictor of first  step = 0  (correspond to real time point 0)
 +
        # Corrector of first  step = 0  (correspond to real time point 1)
 +
        # Corrector of second step = 1  (correspond to real time point 2)
 +
        # Corrector of third  step = 2  (correspond to real time point 3)
 +
        # ...
 +
 
 +
        if curtime == 0:
 +
            SSSidx = 0
 +
        else:
 +
            SSSidx = curtime - 1
 +
 
 +
        # Open output file (coolifc.out0 for first step etc.)
 +
 
 +
        file_in = open('./coolifc.out{:d}'.format(SSSidx),'r')
 +
 
 +
        # Loop over output file to read power distribution
 +
 
 +
        for line in file_in:
 +
            # Split line to values
 +
 
 +
            strtuple = line.split()
 +
 
 +
            # Store power
 +
 
 +
            P.append(float(strtuple[8]))
 +
 
 +
        file_in.close()
 +
 
 +
        if iterating == 0:
 +
            # Moving to next time step
 +
 
 +
            ########################
 +
            # Do some archiving... #
 +
            ########################
 +
 
 +
            # Copy the interface and interface output to
 +
            # archive files for later examination
 +
            # Note: These are Beginning Of Step fields
 +
 
 +
            os.system('cp ./fuel.ifc ./fuel.ifc{:d}'.format(curtime))
 +
 
 +
            # Check if simulation has finished and break out of iterating
 +
            # loop
 +
 
 +
            if (simulating == 0):
 +
                break
 +
 
 +
            ########################################################
 +
            # Update zone burnups (needed for material properties) #
 +
            ########################################################
 +
 
 +
            # Calculate fuel mass of each node (in kg)
 +
 
 +
            m = (math.pi*0.4335**2*20.0)*10.424*1e-3
 +
 
 +
            # Calculate initial HM mass (in kg)
 +
            # (multiply fuel mass by U wt. fraction)
 +
 
 +
            mHM = m*(0.86563+0.015867)
 +
 
 +
            # Get length of current time step (days)
 +
 
 +
            dt = steplengths[curtime]
 +
 
 +
            # Calculate burnup increment for each node
 +
 
 +
            dBU = []
 +
 
 +
            for power in P:
 +
                # Assume constant power throughout the BU step
 +
 
 +
                dBU.append( ((power/1e6)*dt) / mHM) # in MWd/kgU
 +
 
 +
            # Increment burnup of each node
 +
 
 +
            for i in range(10):
 +
                BU[i] = BU[i] + dBU[i]
 +
 
 +
            #######################
 +
            # Increment time step #
 +
            #######################
 +
 
 +
            curdays += steplengths[curtime]
 +
            curtime += 1
 +
 
 +
 
 +
 
 +
        #########################################
 +
        # Calculate fuel temperature solution  #
 +
        #########################################
 +
 
 +
        print("\nMSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS")
 +
        print("|                                          |")
 +
        print("|  Solving fuel temperature for t = {:1d} d  |".format(curdays))
 +
        print("|                                          |")
 +
        print("MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS")
 +
 
 +
        for i in range(10):
 +
 
 +
            ################################################################
 +
            # A polynomial expansion for fuel temperature as a function of #
 +
            # linear power and burnup                                      #
 +
            # and burnup                                                  #
 +
            # Just a stupid power dependence so we don't need to write a  #
 +
            # proper solver for fuel behavior                              #
 +
            ################################################################
 +
 
 +
            # T(bu, P) = a(bu)*LHR + b(bu)
 +
 
 +
            buVec = [0.0,    10.0]
 +
            aVec  = [2.0,    4.0]
 +
            bVec  = [550.0, 550.0]
 +
 
 +
          # Interpolate polynomial coefficients wrt. burnup
 +
 
 +
            a = aVec[0] + (BU[i] - buVec[0])/(buVec[1] - buVec[0])*(aVec[1] -
 +
                                                                    aVec[0])
 +
            b = bVec[0] + (BU[i] - buVec[0])/(buVec[1] - buVec[0])*(bVec[1] -
 +
                                                                    bVec[0])
 +
 
 +
            # Calculate temperature based on power
 +
 
 +
            TFU[i] = a*(P[i]/20.0) + b
 +
 
 +
        ###########################
 +
        # Update interface        #
 +
        ###########################
 +
 
 +
        file_out = open('./fuel.ifc','w')
 +
 
 +
        # Write the header line (TYPE MAT OUT)
 +
 
 +
        file_out.write('2 fuel 0\n')
 +
 
 +
        # Write the mesh type
 +
 
 +
        file_out.write('1\n')
 +
 
 +
        # Write the mesh size (NX XMIN XMAX NY YMIN YMAX NZ ZMIN ZMAX)
 +
 
 +
        file_out.write('1 -0.75 0.75 1 -0.75 0.75 10 -100 100\n')
 +
 
 +
        # Write updated fuel temperatures
 +
 
 +
        for i in range(10):
 +
            # Use the base density throughout the simulation
 +
            # Write density and temperature at this layer
 +
 
 +
            file_out.write('-10.424 {}\n'.format(TFU[i]))
 +
 
 +
        file_out.close()
 +
 
 +
        ############################
 +
        # Signal Serpent (SIGUSR1) #
 +
        ############################
 +
 
 +
        file_out = open('./com.in','w')
 +
 
 +
        file_out.write(str(signal.SIGUSR1.value))
 +
 
 +
        file_out.close()
 +
 
 +
    ####################################
 +
    # Check if simulation has finished #
 +
    ####################################
 +
 
 +
    if (simulating == 0):
 +
        break
 +
 
 +
    ############################
 +
    # Signal Serpent (SIGUSR1) #
 +
    ############################
 +
 
 +
    file_out = open('./com.in','w')
 +
 
 +
    file_out.write(str(signal.SIGUSR1.value))
 +
 
 +
    file_out.close()</nowiki>
  
 
=== input (Stochastic Implicit Euler) ===
 
=== input (Stochastic Implicit Euler) ===
Line 566: Line 992:
 
%    proceeds otherwise normally. This can be used to check that
 
%    proceeds otherwise normally. This can be used to check that
 
%    the calculation sequence itself runs to completion
 
%    the calculation sequence itself runs to completion
−
set acelib "/home/vvvillehe/XSdata/xsdir0K"
+
 
 
%set declib "<path-here>"
 
%set declib "<path-here>"
 
%set nfylib "<path-here>"
 
%set nfylib "<path-here>"
Line 601: Line 1,027:
 
% --- Axial depletion zone division for the fuel material
 
% --- Axial depletion zone division for the fuel material
  
−
div fuel subz 10 -200.0 200.0
+
div fuel subz 10 -100.0 100.0
  
 
% --- Cladding material:
 
% --- Cladding material:
Line 682: Line 1,108:
 
fuel  8  11.808
 
fuel  8  11.808
 
fuel  9  11.808
 
fuel  9  11.808
−
fuel 10  11.808</nowiki>
+
fuel 10  11.808
 +
 
 +
% --- Disable group constant generation
 +
 
 +
set gcu -1</nowiki>
  
 
== Setup ==
 
== Setup ==
  
−
MSCS has been tested with Python 2.7.12 and Serpent 2.1.27.
+
MSCS has been tested with Python 3.8.13 and Serpent 2.2.0.
  
 
#Save the contents of the MSCS.py file (above) to a file of the same name.
 
#Save the contents of the MSCS.py file (above) to a file of the same name.
Line 701: Line 1,131:
  
 
The output of the Serpent run will be printed to the terminal.
 
The output of the Serpent run will be printed to the terminal.
 +
 +
[[Category:Input]]
 +
[[Category:Tutorials]]

Latest revision as of 21:47, 19 August 2022

The Minimal Serpent Coupling Script (MSCS) for burnup calculations is a short (~ 400 lines with comments) Python program intended to give a minimal working example of a wrapper program that can communicate with Serpent in the coupled burnup calculation mode. MSCS provides a working example of externally coupled multi-physics simulations with Serpent and may be a good starting point for the users that are interested in running such simulations with Serpent.

Description

The coupling script communicates with Serpent using the file based communication mode.

Here, the coupling script calculates the fuel temperature solution itself using a not-very-physical functional expansion for fuel temperature as a function of power and local burnup. In many cases that part of the coupling script should be replaced by writing the input for an external solver, running the external solver and reading the results from the external solver output.

The temperature treatment of the interaction physics is done on-the-fly using the TMS temperature treatment technique for the base cross sections and interpolation of thermal scattering data for the thermal scattering libraries.

The problem solved here is a 200 cm long fuel rod in infinite lattice with axially black boundary conditions. The fuel rod is divided axially into 10 depletion zones, for which the fuel temperature is solved by MSCS.

One of the provided inputs uses the explicit Euler discretization for the depletion solution with stochastic approximation based relaxation applied to the tallied transmutation cross sections and one-group fluxes. The second input uses the Stochastic Implicit Euler's method[1].

Files for explicit Euler

MSCS.py (explicit Euler)

#############################################################
#                                                           #
#          Minimal Serpent Coupling Script v 0.3            #
#          For burnup example                               #
#                                                           #
# Created by:  Ville Valtavirta              2016/01/04     #
# Modified by: Ville Valtavirta              2018/10/26     #
#                                                           #
#############################################################

import os
import signal
import math
import time

# Path to Serpent executable

sssexe = '/home/vvvillehe/Serpent2/2.1.31/sss2'

#######################################################
# Create the Serpent input-file for this run          #
# (process id or communication file must be appended) #
#######################################################

# Open original input for reading

file_in = open('./input','r')

# Open a new input file for writing

file_out = open('./coupledinput','w')

# Write original input to new file

for line in file_in:
    file_out.write(line)

# Close original input file

file_in.close()

# Append signalling mode

file_out.write('\n')
file_out.write('set comfile com.in com.out\n')

# Append interface names

file_out.write('\n')
file_out.write('ifc cool.ifc\n\n')
file_out.write('ifc fuel.ifc\n')

# Close new input file

file_out.close()

##############################################
# Write the initial coolant interface file   #
# (Coolant conditions will be held constant) #
##############################################

file_out = open('./cool.ifc','w')

# Write the header line (TYPE MAT OUT)

file_out.write('2 cool 1\n')

# Write the output line (OUTFILE NZ ZMIN ZMAX NR)

file_out.write('coolifc.out 10 -100 100 1\n')

# Write the mesh type

file_out.write('1\n')

# Write the mesh size (NX XMIN XMAX NY YMIN YMAX NZ ZMIN ZMAX)

file_out.write('1 -0.75 0.75 1 -0.75 0.75 1 -100 100\n')

# Write initial coolant temperatures and densities

for i in range(1):
    file_out.write('-0.9  520.0\n')

# Close interface file

file_out.close()

##############################################
# Write the initial fuel interface file      #
# (Fuel temperature will be updated)         #
##############################################

file_out = open('./fuel.ifc','w')

# Write the header line (TYPE MAT OUT)

file_out.write('2 fuel 0\n')

# Write the mesh type

file_out.write('1\n')

# Write the mesh size (NX XMIN XMAX NY YMIN YMAX NZ ZMIN ZMAX)

file_out.write('1 -0.75 0.75 1 -0.75 0.75 10 -100 100\n')

# Write initial fuel temperatures and densities

for i in range(10):
    file_out.write('-10.424 310.0\n')

# Close interface file

file_out.close()

# Archive the initial fuel interface

os.system('cp ./fuel.ifc ./fuel.ifc0')

################################
# Start the Serpent simulation #
################################

# Create a command string that will start the Serpent simulation

runcommand = sssexe+' -omp 3 ./coupledinput &'

# Execute the command string

os.system(runcommand)

#########################################
# Initialize the fuel behavior solution #
#########################################

TFU = []
BU  = []

for i in range(10):
    TFU.append(300.0)
    BU.append(0.0)

# Reset time step

curtime = 0
curdays = 0

# Depletion step lengths (in days)

steplengths = [2, 3]

#############################
# Loop over depletion steps #
#############################

simulating = 1

while simulating == 1:

    ####################################################
    # Picard iteration loop for current depletion step #
    ####################################################

    iterating = 1

    while iterating == 1:
        ###################
        # Wait for signal #
        ###################

        sleeping = 1

        while sleeping == 1:

            # Sleep for two seconds

            time.sleep(2)

            # Open file to check if we got a signal

            fin = open('./com.out','r')

            # Read line

            line = fin.readline()

            # Close file

            fin.close()

            # Check signal

            if int(line) != -1:
                if int(line) == signal.SIGUSR1.value:
                    # Got the signal to resume

                    sleeping = 0

                elif int(line) == signal.SIGUSR2.value:
                    # Got the signal to move to next time point

                    iterating = 0
                    sleeping = 0
                elif int(line) == signal.SIGTERM.value:
                    # Got the signal to end the calculation

                    iterating = 0
                    sleeping = 0
                    simulating = 0
                else:
                    # Unknown signal

                    print("\nUnknown signal read from file, exiting\n")

                    # Exit

                    quit()

                # Reset the signal in the file

                file_out = open('./com.out','w')

                file_out.write('-1')

                file_out.close()

        ###########################
        # Read power distribution #
        ###########################

        # Reset power distribution

        P = []

        # Open output file (coolifc.out0 for first step etc.)

        file_in = open('./coolifc.out{:d}'.format(curtime),'r')

        # Loop over output file to read power distribution

        for line in file_in:
            # Split line to values

            strtuple = line.split()

            # Store power

            P.append(float(strtuple[8]))

        file_in.close()

        if iterating == 0:
            # Moving to next time step

            ########################
            # Do some archiving... #
            ########################

            # Copy the interface and interface output to
            # archive files for later examination
            # Note: These are Beginning Of Step fields

            os.system('cp ./fuel.ifc ./fuel.ifc{:d}'.format(curtime))

            # Check if simulation has finished and break out of iterating
            # loop

            if (simulating == 0):
                break

            ########################################################
            # Update zone burnups (needed for material properties) #
            ########################################################

            # Calculate fuel mass of each node (in kg)

            m = (math.pi*0.4335**2*20.0)*10.424*1e-3

            # Calculate initial HM mass (in kg)
            # (multiply fuel mass by U wt. fraction)

            mHM = m*(0.86563+0.015867)

            # Get length of current time step (days)

            dt = steplengths[curtime]

            # Calculate burnup increment for each node

            dBU = []

            for power in P:
                # Assume constant power throughout the BU step

                dBU.append( ((power/1e6)*dt) / mHM) # in MWd/kgU

            # Increment burnup of each node

            for i in range(10):
                BU[i] = BU[i] + dBU[i]

            #######################
            # Increment time step #
            #######################

            curdays += steplengths[curtime]
            curtime += 1



        #########################################
        # Calculate fuel temperature solution   #
        #########################################

        print("\nMSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS")
        print("|                                          |")
        print("|   Solving fuel temperature for t = {:1d} d   |".format(curdays))
        print("|                                          |")
        print("MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS")

        for i in range(10):

            ################################################################
            # A polynomial expansion for fuel temperature as a function of #
            # linear power and burnup                                      #
            # and burnup                                                   #
            # Just a stupid power dependence so we don't need to write a   #
            # proper solver for fuel behavior                              #
            ################################################################

            # T(bu, P) = a(bu)*LHR + b(bu)

            buVec = [0.0,    10.0]
            aVec  = [2.0,     4.0]
            bVec  = [550.0, 550.0]

           # Interpolate polynomial coefficients wrt. burnup

            a = aVec[0] + (BU[i] - buVec[0])/(buVec[1] - buVec[0])*(aVec[1] -
                                                                    aVec[0])
            b = bVec[0] + (BU[i] - buVec[0])/(buVec[1] - buVec[0])*(bVec[1] -
                                                                    bVec[0])

            # Calculate temperature based on power

            TFU[i] = a*(P[i]/20.0) + b

        ###########################
        # Update interface        #
        ###########################

        file_out = open('./fuel.ifc','w')

        # Write the header line (TYPE MAT OUT)

        file_out.write('2 fuel 0\n')

        # Write the mesh type

        file_out.write('1\n')

        # Write the mesh size (NX XMIN XMAX NY YMIN YMAX NZ ZMIN ZMAX)

        file_out.write('1 -0.75 0.75 1 -0.75 0.75 10 -100 100\n')

        # Write updated fuel temperatures

        for i in range(10):
            # Use the base density throughout the simulation
            # Write density and temperature at this layer

            file_out.write('-10.424 {}\n'.format(TFU[i]))

        file_out.close()

        ############################
        # Signal Serpent (SIGUSR1) #
        ############################

        file_out = open('./com.in','w')

        file_out.write(str(signal.SIGUSR1.value))

        file_out.close()

    ####################################
    # Check if simulation has finished #
    ####################################

    if (simulating == 0):
        break

    ############################
    # Signal Serpent (SIGUSR1) #
    ############################

    file_out = open('./com.in','w')

    file_out.write(str(signal.SIGUSR1.value))

    file_out.close()

input (explicit Euler)

% --- Input for MSCS testing

set title "Serpent-MSCS externally coupled burnup calculation"

%set acelib "<path-here>"

% --- Tip: In order to quickly obtain unphysical results you can
%     comment out the declib and nfylib. This means that fission products
%     and actinides will not be produced, but the calculation sequence
%     proceeds otherwise normally. This can be used to check that
%     the calculation sequence itself runs to completion

%set declib "<path-here>"
%set nfylib "<path-here>"

% --- Fuel Pin definition:

pin 1
fuel   0.4335
gas    0.442000
clad   0.502500
cool

% --- Lattice (type = 1, pin pitch = 1.5):

lat 10  1  0.0 0.0 1 1 1.5
1

% --- Boundary of the geometry (200 cm active length)

surf 2 cuboid -0.75 0.75 -0.75 0.75 -100 100

% --- Cell definitions:

cell  3  0  fill 10  -2     % Pin-cell
cell 99  0  outside   2     % Outside world

% --- Fuel material:

mat fuel    -10.424 tft 309 3000.0 burn 1 rgb 100 160 140
 92235.03c   -0.015867
 92238.03c   -0.86563
  8016.03c   -0.1185

% --- Axial depletion zone division for the fuel material

div fuel subz 10 -100.0 100.0

% --- Cladding material:

mat clad     -6.55 rgb 200 200 200
 40090.03c   -0.98135
 24052.03c   -0.00100
 26056.03c   -0.00135
 28058.03c   -0.00055
 50120.03c   -0.01450
  8016.03c   -0.00125

% --- Gas gap:

mat gas -1E-4 rgb 255 255 255
 2004.06c     1.0

% --- Coolant (Temperature won't change, but will be given via interface):

mat cool     -0.90 moder lwtr 1001 rgb 150 160 240
 1001.03c     0.66667
 8016.03c     0.33333
 5010.03c     0.0001

% --- Thermal scattering data for light water:
%     On-the-fly treatment for SAB-data between 474 K -- 624 K

therm lwtr 520 lwj3.07t lwj3.09t

% --- Axially black boundary condition

set bc 3 3 1

% --- System power corresponding to a mean linear power of 200 W/cm

set power 40000.0

% --- Simulation population (very small -> non-physical results)

set pop 5000 50 50

% --- Maximum number of coupled calculation iterations set to 2

set ccmaxiter 2

% --- Depletion scheme constant extrapolation == explicit Euler

set pcc ce

% --- Depletion history (just two steps)

dep daystep 2 3

% --- Geometry plot (xz)

plot 2 500 1000

% --- xy-plot of fission power / thermal flux

mesh 3 500 500

% --- xz-plot of fission power / thermal flux

mesh 2 500 1000

% --- xy-plot of temperature distribution / thermal flux

mesh 10 3 500 500

% --- xz-plot of temperature distribution / thermal flux

mesh 10 2 500 1000

% --- Set material volumes for depletion zones

set mvol
fuel  1   11.808
fuel  2   11.808
fuel  3   11.808
fuel  4   11.808
fuel  5   11.808
fuel  6   11.808
fuel  7   11.808
fuel  8   11.808
fuel  9   11.808
fuel 10   11.808

% --- Disable group constant generation

set gcu -1

Files for stochastic implicit Euler

MSCS.py (Stochastic Implicit Euler)

#############################################################
#                                                           #
#          Minimal Serpent Coupling Script v 0.2            #
#          For burnup example (Stochastic Implicit Euler)   #
#                                                           #
# Created by:  Ville Valtavirta              2016/01/04     #
# Modified by: Ville Valtavirta              2020/10/26     #
#                                                           #
#############################################################

import os
import signal
import math
import time

# Path to Serpent executable

sssexe = '/home/vvvillehe/Serpent2/2.1.31/sss2'

#######################################################
# Create the Serpent input-file for this run          #
# (process id or communication file must be appended) #
#######################################################

# Open original input for reading

file_in = open('./input','r')

# Open a new input file for writing

file_out = open('./coupledinput','w')

# Write original input to new file

for line in file_in:
    file_out.write(line)

# Close original input file

file_in.close()

# Append signalling mode

file_out.write('\n')
file_out.write('set comfile com.in com.out\n')

# Append interface names

file_out.write('\n')
file_out.write('ifc cool.ifc\n\n')
file_out.write('ifc fuel.ifc\n')

# Close new input file

file_out.close()

##############################################
# Write the initial coolant interface file   #
# (Coolant conditions will be held constant) #
##############################################

file_out = open('./cool.ifc','w')

# Write the header line (TYPE MAT OUT)

file_out.write('2 cool 1\n')

# Write the output line (OUTFILE NZ ZMIN ZMAX NR)

file_out.write('coolifc.out 10 -100 100 1\n')

# Write the mesh type

file_out.write('1\n')

# Write the mesh size (NX XMIN XMAX NY YMIN YMAX NZ ZMIN ZMAX)

file_out.write('1 -0.75 0.75 1 -0.75 0.75 1 -100 100\n')

# Write initial coolant temperatures and densities

for i in range(1):
    file_out.write('-0.9  520.0\n')

# Close interface file

file_out.close()

##############################################
# Write the initial fuel interface file      #
# (Fuel temperature will be updated)         #
##############################################

file_out = open('./fuel.ifc','w')

# Write the header line (TYPE MAT OUT)

file_out.write('2 fuel 0\n')

# Write the mesh type

file_out.write('1\n')

# Write the mesh size (NX XMIN XMAX NY YMIN YMAX NZ ZMIN ZMAX)

file_out.write('1 -0.75 0.75 1 -0.75 0.75 10 -100 100\n')

# Write initial fuel temperatures and densities

for i in range(10):
    file_out.write('-10.424 310.0\n')

# Close interface file

file_out.close()

# Archive the initial fuel interface

os.system('cp ./fuel.ifc ./fuel.ifc0')

################################
# Start the Serpent simulation #
################################

# Create a command string that will start the Serpent simulation

runcommand = sssexe+' -omp 3 ./coupledinput &'

# Execute the command string

os.system(runcommand)

#########################################
# Initialize the fuel behavior solution #
#########################################

TFU = []
BU  = []

for i in range(10):
    TFU.append(300.0)
    BU.append(0.0)

# Reset time step

curtime = 0
curdays = 0

# Depletion step lengths (in days)

steplengths = [2, 3]

#############################
# Loop over depletion steps #
#############################

simulating = 1

while simulating == 1:

    ####################################################
    # Picard iteration loop for current depletion step #
    ####################################################

    iterating = 1

    while iterating == 1:
        ###################
        # Wait for signal #
        ###################

        sleeping = 1

        while sleeping == 1:

            # Sleep for two seconds

            time.sleep(2)

            # Open file to check if we got a signal

            fin = open('./com.out','r')

            # Read line

            line = fin.readline()

            # Close file

            fin.close()

            # Check signal

            if int(line) != -1:
                if int(line) == signal.SIGUSR1.value:
                    # Got the signal to resume

                    sleeping = 0

                elif int(line) == signal.SIGUSR2.value:
                    # Got the signal to move to next time point

                    iterating = 0
                    sleeping = 0
                elif int(line) == signal.SIGTERM.value:
                    # Got the signal to end the calculation

                    iterating = 0
                    sleeping = 0
                    simulating = 0
                else:
                    # Unknown signal

                    print("\nUnknown signal read from file, exiting\n")

                    # Exit

                    quit()

                # Reset the signal in the file

                file_out = open('./com.out','w')

                file_out.write('-1')

                file_out.close()

        ###########################
        # Read power distribution #
        ###########################

        # Reset power distribution

        P = []

        # Get output file number corresponding for current time
        # Serpent versions up to 2.1.29 use the following numbering with SIE:
        # Predictor of first  step = 0  (correspond to real time point 0)
        # Corrector of first  step = 0  (correspond to real time point 1)
        # Corrector of second step = 1  (correspond to real time point 2)
        # Corrector of third  step = 2  (correspond to real time point 3)
        # ...

        if curtime == 0:
            SSSidx = 0
        else:
            SSSidx = curtime - 1

        # Open output file (coolifc.out0 for first step etc.)

        file_in = open('./coolifc.out{:d}'.format(SSSidx),'r')

        # Loop over output file to read power distribution

        for line in file_in:
            # Split line to values

            strtuple = line.split()

            # Store power

            P.append(float(strtuple[8]))

        file_in.close()

        if iterating == 0:
            # Moving to next time step

            ########################
            # Do some archiving... #
            ########################

            # Copy the interface and interface output to
            # archive files for later examination
            # Note: These are Beginning Of Step fields

            os.system('cp ./fuel.ifc ./fuel.ifc{:d}'.format(curtime))

            # Check if simulation has finished and break out of iterating
            # loop

            if (simulating == 0):
                break

            ########################################################
            # Update zone burnups (needed for material properties) #
            ########################################################

            # Calculate fuel mass of each node (in kg)

            m = (math.pi*0.4335**2*20.0)*10.424*1e-3

            # Calculate initial HM mass (in kg)
            # (multiply fuel mass by U wt. fraction)

            mHM = m*(0.86563+0.015867)

            # Get length of current time step (days)

            dt = steplengths[curtime]

            # Calculate burnup increment for each node

            dBU = []

            for power in P:
                # Assume constant power throughout the BU step

                dBU.append( ((power/1e6)*dt) / mHM) # in MWd/kgU

            # Increment burnup of each node

            for i in range(10):
                BU[i] = BU[i] + dBU[i]

            #######################
            # Increment time step #
            #######################

            curdays += steplengths[curtime]
            curtime += 1



        #########################################
        # Calculate fuel temperature solution   #
        #########################################

        print("\nMSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS")
        print("|                                          |")
        print("|   Solving fuel temperature for t = {:1d} d   |".format(curdays))
        print("|                                          |")
        print("MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS-MSCS")

        for i in range(10):

            ################################################################
            # A polynomial expansion for fuel temperature as a function of #
            # linear power and burnup                                      #
            # and burnup                                                   #
            # Just a stupid power dependence so we don't need to write a   #
            # proper solver for fuel behavior                              #
            ################################################################

            # T(bu, P) = a(bu)*LHR + b(bu)

            buVec = [0.0,    10.0]
            aVec  = [2.0,     4.0]
            bVec  = [550.0, 550.0]

           # Interpolate polynomial coefficients wrt. burnup

            a = aVec[0] + (BU[i] - buVec[0])/(buVec[1] - buVec[0])*(aVec[1] -
                                                                    aVec[0])
            b = bVec[0] + (BU[i] - buVec[0])/(buVec[1] - buVec[0])*(bVec[1] -
                                                                    bVec[0])

            # Calculate temperature based on power

            TFU[i] = a*(P[i]/20.0) + b

        ###########################
        # Update interface        #
        ###########################

        file_out = open('./fuel.ifc','w')

        # Write the header line (TYPE MAT OUT)

        file_out.write('2 fuel 0\n')

        # Write the mesh type

        file_out.write('1\n')

        # Write the mesh size (NX XMIN XMAX NY YMIN YMAX NZ ZMIN ZMAX)

        file_out.write('1 -0.75 0.75 1 -0.75 0.75 10 -100 100\n')

        # Write updated fuel temperatures

        for i in range(10):
            # Use the base density throughout the simulation
            # Write density and temperature at this layer

            file_out.write('-10.424 {}\n'.format(TFU[i]))

        file_out.close()

        ############################
        # Signal Serpent (SIGUSR1) #
        ############################

        file_out = open('./com.in','w')

        file_out.write(str(signal.SIGUSR1.value))

        file_out.close()

    ####################################
    # Check if simulation has finished #
    ####################################

    if (simulating == 0):
        break

    ############################
    # Signal Serpent (SIGUSR1) #
    ############################

    file_out = open('./com.in','w')

    file_out.write(str(signal.SIGUSR1.value))

    file_out.close()

input (Stochastic Implicit Euler)

% --- Input for MSCS testing

set title "Serpent-MSCS externally coupled burnup calculation"

%set acelib "<path-here>"

% --- Tip: In order to quickly obtain unphysical results you can
%     comment out the declib and nfylib. This means that fission products
%     and actinides will not be produced, but the calculation sequence
%     proceeds otherwise normally. This can be used to check that
%     the calculation sequence itself runs to completion

%set declib "<path-here>"
%set nfylib "<path-here>"

% --- Fuel Pin definition:

pin 1
fuel   0.4335
gas    0.442000
clad   0.502500
cool

% --- Lattice (type = 1, pin pitch = 1.5):

lat 10  1  0.0 0.0 1 1 1.5
1

% --- Boundary of the geometry (200 cm active length)

surf 2 cuboid -0.75 0.75 -0.75 0.75 -100 100

% --- Cell definitions:

cell  3  0  fill 10  -2     % Pin-cell
cell 99  0  outside   2     % Outside world

% --- Fuel material:

mat fuel    -10.424 tft 309 3000.0 burn 1 rgb 100 160 140
 92235.03c   -0.015867
 92238.03c   -0.86563
  8016.03c   -0.1185

% --- Axial depletion zone division for the fuel material

div fuel subz 10 -100.0 100.0

% --- Cladding material:

mat clad     -6.55 rgb 200 200 200
 40090.03c   -0.98135
 24052.03c   -0.00100
 26056.03c   -0.00135
 28058.03c   -0.00055
 50120.03c   -0.01450
  8016.03c   -0.00125

% --- Gas gap:

mat gas -1E-4 rgb 255 255 255
 2004.06c     1.0

% --- Coolant (Temperature won't change, but will be given via interface):

mat cool     -0.90 moder lwtr 1001 rgb 150 160 240
 1001.03c     0.66667
 8016.03c     0.33333
 5010.03c     0.0001

% --- Thermal scattering data for light water:
%     On-the-fly treatment for SAB-data between 474 K -- 624 K

therm lwtr 520 lwj3.07t lwj3.09t

% --- Axially black boundary condition

set bc 3 3 1

% --- System power corresponding to a mean linear power of 200 W/cm

set power 40000.0

% --- Simulation population (very small -> non-physical results)

set pop 5000 50 50

% --- Depletion scheme stochastic implicit Euler with 2 iterations per step

set sie 2

% --- Depletion history (just two steps)

dep daystep 2 3

% --- Geometry plot (xz)

plot 2 500 1000

% --- xy-plot of fission power / thermal flux

mesh 3 500 500

% --- xz-plot of fission power / thermal flux

mesh 2 500 1000

% --- xy-plot of temperature distribution / thermal flux

mesh 10 3 500 500

% --- xz-plot of temperature distribution / thermal flux

mesh 10 2 500 1000

% --- Set material volumes for depletion zones

set mvol
fuel  1   11.808
fuel  2   11.808
fuel  3   11.808
fuel  4   11.808
fuel  5   11.808
fuel  6   11.808
fuel  7   11.808
fuel  8   11.808
fuel  9   11.808
fuel 10   11.808

% --- Disable group constant generation

set gcu -1

Setup

MSCS has been tested with Python 3.8.13 and Serpent 2.2.0.

  1. Save the contents of the MSCS.py file (above) to a file of the same name.
  2. Replace the absolute path to the Serpent executable on line 18 of MSCS.py.
  3. Save the contents of one of the input files to a file called input in the same directory where the MSCS.py file is located.
  4. Add the path for the cross section libraries to the input-file using the set acelib input option.
  5. Add the path for fission yield data to the input-file using the set nfylib input option.
  6. Add the path for radioactive decay data to the input-file using the set declib input option.
  7. The test case is now ready to run.

Note: Before running the MSCS you should familiarize yourself with finding and killing background processes in your operating system. Since MSCS runs Serpent as a background process, killing MSCS will not kill the Serpent process automatically.

Run the simulation from the folder containing both files using the command python MSCS.py

The output of the Serpent run will be printed to the terminal.
  1. ^ J. Dufek and H. Anglart, "Derivation of a stable coupling scheme for Monte Carlo burnup calculations with the thermal-hydraulic feedback", Ann. Nucl. Energy 62, pp. 260 (2013) DOI: 10.1016/j.anucene.2013.06.025
Retrieved from "https://serpent.vtt.fi/mediawiki/index.php?title=MSCS_for_burnup_calculations&oldid=5754"