Reformulate web presence

.github/workflows/testcommands.yml

@@ -0,0 +1,42 @@
+name: test commands
+on:
+  workflow_dispatch:
+  schedule:
+    - cron: '44 4 * * 0'
+
+jobs:
+  test_commands:
+    strategy:
+      matrix:
+        command: ['k_eff','energy']
+    runs-on: ubuntu-22.04
+    steps:
+      - name: checkout_main_repo
+        uses: actions/checkout@v5
+        with:
+          repository: 'openmsr/are'
+          ref: '${{github.sha}}'
+          path: are
+
+      - name: bootstrap
+        run: sudo apt-get -y update
+
+      - name: environment
+        uses: mamba-org/setup-micromamba@v1
+        with:
+          environment-name: cad2omc
+          create-args: >-
+            python=3.10
+            openmc_0.13.3=dagmc_mpi_openmpi*
+          init_shell: bash
+
+      - name: check openmc
+        run: |
+          openmc --version
+          python -c "import openmc"
+      - name: install cad-to-openmc
+        run: |
+          python -m pip install upgrade pip
+          python -m pip install cad_to_openmc
+      - name: run ${{ command }}
+        run: python run.py ${{ command }}

.gitignore

@@ -127,3 +127,12 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+# data-files
+*.dat
+*.txt
+*.h5
+*.png
+*.out
+*.err
+h5m_files/*

README.md

@@ -5,7 +5,7 @@ detailed cad model and simulations of the [are](https://en.wikipedia.org/wiki/Ai
 
 ## model
 
-### are core 
+### are core
 ![](figures/core.png)
 
 [docs/are.pdf](docs/are.pdf) lists reference of the are core design, documented in the ornl reports and located in the repository [github.com/openmsr/msr-archive](https://github.com/openmsr/msr-archive/blob/master/README.md)
@@ -14,24 +14,86 @@ work-in-progress cad model of the are can be found [here](https://cad.onshape.co
 
 note that this work and the cad model is under the GNU General Public License v3.0
 
-## prerequisites 
+## prerequisites
 ### CAD_to_openMC
 [CAD_to_openMC](https://github.com/openmsr/CAD_to_openMC) is an open-source package to convert CAD geometry (in the form of '.step' files) into an openmc-readable h5m file
 
+
+
 ### openmc
-these simulations use [openmc](https://docs.openmc.org/en/stable/). automated source installation scripts for linux can be found [here](https://github.com/openmsr/openmc_install_scripts)
+These simulations use [openmc](https://docs.openmc.org/en/stable/). automated source installation scripts for linux can be found [here](https://github.com/openmsr/openmc_install_scripts)
+
+
 
 ## simulation guide
+We have set up a set of helper scripts to run a set of standard simulations
+with the ARE-geometry: k_eff, neutron_flux, ...
+
+All these runs first check if the necessary set of three .h5m-files exist. If
+they do, they are used, if not a set will be created on the fly using the helper
+script are_step_to_h5m.py.  Thus, should you wish to regenerate the set (or
+part of it) you can simply move the files out of the way or run
+are_step_to_h5m.py manually.
+
+### k\_eff
+This is often the first simulation run done on a geometry which
+estimates the reactivity or criticality of the reactor geometry.
+
+### energy
+This run defines an overall tally in the reactor geometry that measures the
+energy of neutron within it, as well as a set of localized energy-scoring
+tallies only scoring particles within:
+1. fuel-salt within the reactor core
+2. the coolant salt.
+
+Please note that even though the ARE reactor geometry describes a fluid
+fuel-reactor, this simulation is a static steady-state run, meaning we
+presently neglect depletion/burn-up and things like delayed neutrons.
+It is possible to account for such effects, but is considere out-of-scope for
+this set of sims.
+
+### neutron_flux
+This simulation run is similar to k_eff, ecept that it adds a set of
+mesh-tallies on slices through the reactor core. Each plane is set
+to score the neutron-flux in the regions as a function of spatial coordinate,
+or more precisely flux binnned in cuboid volume elements in the planes.
+The tallies predefined in the computations are a set of
+XY|z\in{-50,-25,0,25,50} slices as well as slices at XZ|y=0, and YZ|z=0.
+Tallies, and mesh-tallies in particular, create some extra overhead slowing the
+simulations down somewhat - this is normal.
+
+### photon_flux
+Very similar to the neutron flux run except in this case what is monitored is the
+photon flux, i.e. the gamma-radiation also produced in the reactor.
+
+### geometry plot
+This type of run is not a simulation run as such, instead it creates a set of
+geometry slice plots aælong the principal axis. This is useful to make sure
+that you are in fact simulating the geometry you think you are.
+Lately, a new tool has been included the openmc-distribution: openmc-plotter,
+which largely renders doing this type of run manually moot. The plotter does perform a geometry plot-run in the background.
+
+### volume
+Run a stochastic volume calculation on the ARE-geometry. This type of computation estimates the volume occupied by each defined
+material in the model.
 
-first, clone the repository
+### crit_approach
+One of the 1st experiments performed with the ARE was to approach a the
+cirticality point of the reactor. This was done by adding small batches of
+fissile fuel to the carrier salt in the reactor. Using the experimental logs of
+the time (see ORNL-) we have estimated the amount of fuel added and hence the
+constituent concentration of various nuclides in the ARE-core.
+The procedure was to gradually increase the amount of fissile material all the while
+measuring reactivity, until criticality was reached.
+Ths simulation sript mimics this procedure and plots the result for 35 in. withdrawn.
 
-```
-git clone https://github.com/openmsr/are.git
-```
+ref: Cottrell, W. B., et. al, "Operation of the Aircraft Reactor Experiment", ORNL-1845, 1955
 
-enter the are folder and run the `run.sh` script
+### shim_rod_calibration
+In the latter stages of the criticality approach experiment, further measurements were
+performed; The experimental procedure was to insert the safety rods half-way (20 in. in), add fuel, and
+then measure the reactivity at set of safety rod positions: 20,25,30,35 in. withdrawn.
+This way experimenters obtained a first estimate of the safety rod worth.
 
-```
-cd are
-bash run.sh
-```
+This simulation runs attempts to mimic this set of measurements.
+ref: Cottrell, W. B., et. al, "Operation of the Aircraft Reactor Experiment", ORNL-1845, 1955

scripts/are_criticality_approach.py

@@ -0,0 +1,57 @@
+import os
+from fuel_salts import *
+from build_model import ARE_reactor
+import matplotlib.pyplot as plt
+
+###############################################################################
+# ARE criticality approach experiment
+###############################################################################
+
+# inputs (nuclear data path is optional)
+h5m_files=['h5m_files/ARE_core.h5m','h5m_files/ARE_safety_rod.h5m','h5m_files/ARE_regulating_rod.h5m']
+print(h5m_files)
+
+shim_rod_z = 35.0
+ctrl_rod_z = 0.0
+run_labels = ( '1','2','3','4','5','6','7','8d','9b','10c','11c','12b' )
+crit_run_salts = [ salts[l] for l in run_labels ]
+
+# writes k_eff & stdev separated by a space on each line of the output file
+output_filename = 'k_effs_criticality_approach.txt'
+with open(output_filename,'w') as k_file:
+    print(f'#shim_rod_z={shim_rod_z} in, ctrl_rod_z={ctrl_rod_z} in', file=k_file)
+    print('#k_eff k_eff_stddev', file=k_file)
+
+for f in crit_run_salts:
+    if f.get_nuclide_densities()['U235'].percent<=0:
+      with open(output_filename,'a') as k_file:
+        k_file.write('0.0 0.0\n')
+    else:
+      reactor=ARE_reactor(h5m_files[0],h5m_files[1],h5m_files[2],f,shim_rod_z=[shim_rod_z,shim_rod_z,shim_rod_z],ctrl_rod_z=ctrl_rod_z,batches=5,inactive=1, particles=1000)
+      #reactor.bld_model()
+      sp_filepath = reactor.model.run(output = True)
+      with openmc.StatePoint(sp_filepath) as sp, open(output_filename,'a') as k_file:
+        k_eff = sp.keff
+        k_file.write(f"{k_eff.nominal_value} {k_eff.std_dev}\n")
+
+# re-read the data-file
+k=np.ndarray([len(crit_run_salts),2])
+with open(output_filename,'r') as f:
+    i=0
+    for line in f:
+        if line[0]=='#':
+            continue
+        k[i,:]=line.split()
+        i=i+1
+print(k)
+
+U235lbs=[concentrations[l][3]/100*concentrations[l][8] for l in run_labels]
+
+f,a = plt.subplots()
+a.errorbar(U235lbs, k[:,0], yerr=k[:,1], capsize=3,label=f'shim_rods at {shim_rod_z} in.')
+a.set_xlabel('U235 inventory [lbs]')
+a.set_ylabel('k_{eff} [ ]')
+a.grid()
+plt.tight_layout()
+plt.legend()
+plt.savefig(f'k_eff_vs_U_sz{shim_rod_z}.png')

scripts/are_energy.py

@@ -0,0 +1,119 @@
+import openmc
+import matplotlib.pyplot as plt
+import numpy as np
+
+from build_model import ARE_reactor
+
+###############################################################################
+#create .png files of photon flux (all safety rods inserted)
+###############################################################################
+
+ARE = ARE_reactor('h5m_files/ARE_core.h5m','h5m_files/ARE_safety_rod.h5m','h5m_files/ARE_regulating_rod.h5m')
+ARE.bld_model()
+
+#tallies
+tallies = openmc.Tallies()
+
+# particle filters
+neutron_particle_filter=openmc.ParticleFilter(['neutron'])
+photon_particle_filter=openmc.ParticleFilter(['photon'])
+
+# energy binned filters
+energy_bins = openmc.mgxs.GROUP_STRUCTURES['VITAMIN-J-175']
+energy_filter = openmc.EnergyFilter(energy_bins)
+
+#domain filters
+print([m.id for m in ARE.materials if m.name=='salt'])
+fuel_domain = openmc.MaterialFilter([m.id for m in ARE.materials if m.name=='salt'])
+coolant_domain = openmc.MaterialFilter([m.id for m in ARE.materials if m.name=='coolant'])
+
+# tallies
+tallyEn0 = openmc.Tally(name='energy_n')
+tallyEn0.filters = [energy_filter,neutron_particle_filter]
+tallyEn0.scores = ['flux']
+tallies.append(tallyEn0)
+
+tallyEp0 = openmc.Tally(name='energy_p')
+tallyEp0.filters = [energy_filter,photon_particle_filter]
+tallyEp0.scores = ['flux']
+tallies.append(tallyEp0)
+
+tallyEn1 = openmc.Tally(name='energy_fuel_n')
+tallyEn1.filters = [energy_filter,neutron_particle_filter,fuel_domain]
+tallyEn1.scores = ['flux']
+tallies.append(tallyEn1)
+
+tallyEp1 = openmc.Tally(name='energy_fuel_p')
+tallyEp1.filters = [energy_filter,photon_particle_filter,fuel_domain]
+tallyEp1.scores = ['flux']
+tallies.append(tallyEp1)
+
+tallyEn1 = openmc.Tally(name='energy_coolant_n')
+tallyEn1.filters = [energy_filter,neutron_particle_filter,coolant_domain]
+tallyEn1.scores = ['flux']
+tallies.append(tallyEn1)
+
+tallyEp1 = openmc.Tally(name='energy_coolant_p')
+tallyEp1.filters = [energy_filter,photon_particle_filter,coolant_domain]
+tallyEp1.scores = ['flux']
+tallies.append(tallyEp1)
+
+# override the default set of tallies
+ARE.model.tallies=tallies
+ARE.model.settings.photon_transport=True
+
+sp_filename = ARE.model.run()
+sp = openmc.StatePoint(sp_filename)
+
+n_tally0 = sp.get_tally(name='energy_n')
+p_tally0 = sp.get_tally(name='energy_p')
+n_tallyf = sp.get_tally(name='energy_fuel_n')
+p_tallyf = sp.get_tally(name='energy_fuel_p')
+n_tallyc = sp.get_tally(name='energy_coolant_n')
+p_tallyc = sp.get_tally(name='energy_coolant_p')
+
+
+fluxn0 = n_tally0.get_slice(scores=['flux'])
+fluxp0 = p_tally0.get_slice(scores=['flux'])
+fluxnf = n_tallyf.get_slice(scores=['flux'])
+fluxpf = p_tallyf.get_slice(scores=['flux'])
+fluxnc = n_tallyc.get_slice(scores=['flux'])
+fluxpc = p_tallyc.get_slice(scores=['flux'])
+
+energy_midpts=(energy_bins[:-1]+energy_bins[1:])*0.5
+
+fig,ax = plt.subplots()
+pos=ax.errorbar(energy_midpts,fluxn0.mean.squeeze(), yerr=fluxn0.std_dev.squeeze(), capsize=3, capthick=1, label='neutron flux')
+pos=ax.errorbar(energy_midpts,fluxp0.mean.squeeze(), yerr=fluxp0.std_dev.squeeze(), capsize=3, capthick=1, label='photon flux')
+ax.grid()
+ax.set_title('flux overall')
+ax.set_xscale('log')
+ax.set_xlabel('energy [eV]')
+ax.set_ylabel('flux [cts. / (s s.p)]')
+plt.tight_layout()
+plt.legend()
+plt.savefig(f'flux_energy', dpi=256)
+
+fig,ax = plt.subplots()
+pos=ax.errorbar(energy_midpts,fluxnf.mean.squeeze(), yerr=fluxnf.std_dev.squeeze(), capsize=3, capthick=1, label='neutron flux')
+pos=ax.errorbar(energy_midpts,fluxpf.mean.squeeze(), yerr=fluxpf.std_dev.squeeze(), capsize=3, capthick=1, label='photon flux')
+ax.grid()
+ax.set_title('flux in fuel salt')
+ax.set_xscale('log')
+ax.set_xlabel('energy [eV]')
+ax.set_ylabel('flux [cts. / (s s.p)]')
+plt.tight_layout()
+plt.legend()
+plt.savefig(f'flux_energy_fuel', dpi=256)
+
+fig,ax = plt.subplots()
+pos=ax.errorbar(energy_midpts,fluxnc.mean.squeeze(), yerr=fluxnc.std_dev.squeeze(), capsize=3, capthick=1, label='neutron flux')
+pos=ax.errorbar(energy_midpts,fluxpc.mean.squeeze(), yerr=fluxpc.std_dev.squeeze(), capsize=3, capthick=1, label='photon flux')
+ax.grid()
+ax.set_title('flux in fuel salt')
+ax.set_xscale('log')
+ax.set_xlabel('energy [eV]')
+ax.set_ylabel('flux [cts. / (s s.p)]')
+plt.tight_layout()
+plt.legend()
+plt.savefig(f'flux_energy_coolant', dpi=256)