Added functionality: reading from an external source file by bonnettheophile · Pull Request #194 · CambridgeNuclear/SCONE

bonnettheophile · 2026-02-09T14:54:27Z

Added class externalSource extending source_inter. Source file can be read in binary or in ASCII, using the formatting of printSource.
Added functionality to printStatement to use binary format instead of ASCII.
In physicsPackages, setting printSource 2; prints to binary format

Modified files:
eigenPhysicsPackage_class.f90
fixedSourcePhysicsPackage_class.f90
sourceFactory_func.f90
source_inter.f90
particleDungeon_class.f90
genericProcedures.f90

Added file: externalSource_class.f90

* Added class externalSource extending source_inter. Source file can be read in binary or in ASCII, using the formatting of printSource. * Added functionality to printStatement to use binary format instead of ASCII. Added a new readBinary switch to physicPackages.

ChasingNeutrons · 2026-02-09T15:15:55Z

A small, quick comment: I would change the name to, perhaps, 'fileSource'. This is because 'external source' is commonly used interchangeably with 'fixed source', so it might be a bit ambiguous.

valeriaRaffuzzi

Thanks a lot for the development, this is very nice!! I just have a few comments and questions. Aside from those I left throughout the code, more high level things:

I do agree with Paul about changing the name to fileSource
We should also update the docs, i.e., the InputManual to include the change.

valeriaRaffuzzi · 2026-02-09T17:18:47Z

ParticleObjects/particleDungeon_class.f90

+    if (writeBinary) then
+      filename = trim(name)//'.bin'
+      open(unit = id, file = filename, status = 'replace', access = 'stream', form = 'unformatted')
+
+      ! Print out each particle co-ordinate
+      do i = 1, self % pop
+        write(id) self % prisoners(i) % r, self % prisoners(i) % dir, &
+                  self % prisoners(i) % E, real(self % prisoners(i) % G, defReal), &
+                  real(self % prisoners(i) % broodID, defReal)
+      end do
+    else
+      ! Print out each particle co-ordinate
+      filename = trim(name)//'.txt'
+      open(unit = id, file = filename, status = 'replace')

-    ! Print out each particle co-ordinate
-    do i = 1, self % pop
-      write(10, *) self % prisoners(i) % r, self % prisoners(i) % dir, &
-                   self % prisoners(i) % E, self % prisoners(i) % G, &
-                   self % prisoners(i) % broodID
-    end do
+      ! Print out each particle co-ordinate
+      do i = 1, self % pop
+        write(id,*) self % prisoners(i) % r, self % prisoners(i) % dir, &
+                    self % prisoners(i) % E, self % prisoners(i) % G, &
+                    self % prisoners(i) % broodID
+      end do
+    end if


Could you have something like this, where the loop over particles is only written once?

Why do you need real() to write to binary and not for ascii?

I think the group should be an integer, since it is an index.

Suggested change

if (writeBinary) then

filename = trim(name)//'.bin'

open(unit = id, file = filename, status = 'replace', access = 'stream', form = 'unformatted')

! Print out each particle co-ordinate

do i = 1, self % pop

write(id) self % prisoners(i) % r, self % prisoners(i) % dir, &

self % prisoners(i) % E, real(self % prisoners(i) % G, defReal), &

real(self % prisoners(i) % broodID, defReal)

end do

else

! Print out each particle co-ordinate

filename = trim(name)//'.txt'

open(unit = id, file = filename, status = 'replace')

! Print out each particle co-ordinate

do i = 1, self % pop

write(10, *) self % prisoners(i) % r, self % prisoners(i) % dir, &

self % prisoners(i) % E, self % prisoners(i) % G, &

self % prisoners(i) % broodID

end do

! Print out each particle co-ordinate

do i = 1, self % pop

write(id,*) self % prisoners(i) % r, self % prisoners(i) % dir, &

self % prisoners(i) % E, self % prisoners(i) % G, &

self % prisoners(i) % broodID

end do

end if

if (writeBinary) then

filename = trim(name)//'.bin'

open(unit = id, file = filename, status = 'replace', access = 'stream', form = 'unformatted')

else

! Print out each particle co-ordinate

filename = trim(name)//'.txt'

open(unit = id, file = filename, status = 'replace')

end if

! Print out each particle co-ordinate

do i = 1, self % pop

write(id) self % prisoners(i) % r, self % prisoners(i) % dir, &

self % prisoners(i) % E, real(self % prisoners(i) % G, defReal), &

real(self % prisoners(i) % broodID, defReal)

end do

valeriaRaffuzzi · 2026-02-09T17:20:45Z

PhysicsPackages/eigenPhysicsPackage_class.f90

      ! Update RNG after it was used to normalise particle population
      call self % pRNG % stride(1)

+      ! Print source in ASCII format if requested


Could also avoid the repeated call by setting a binary flag that's inputted into the function printToFile

valeriaRaffuzzi · 2026-02-09T17:21:27Z

PhysicsPackages/fixedSourcePhysicsPackage_class.f90

      ! Update RNG after source generation
      call self % pRNG % stride(self % totalPop)

+      ! Print source in ASCII if requested


Same comment as for the eigenPP

valeriaRaffuzzi · 2026-02-09T17:22:48Z

SharedModules/genericProcedures.f90

+        output = temp
+      case (iostat_end)
+        EOF = .true.
+        return


Surely you don't need the return, since there is nothing after the select case.

valeriaRaffuzzi · 2026-02-09T17:23:22Z

SharedModules/genericProcedures.f90

  end function printParticleType

+  !!
+  !! Read a line from the source file in ASCII format


The comment shoudl be 'ASCII or binary' I believe?

I would write a comment about what EOF is.

valeriaRaffuzzi · 2026-02-12T16:14:08Z

ParticleObjects/Source/externalSource_class.f90

+    allocate(self % r(3, self % numberNeutrons))
+    allocate(self % dir(3, self % numberNeutrons))


Suggested change

allocate(self % r(3, self % numberNeutrons))

allocate(self % dir(3, self % numberNeutrons))

allocate(self % r(3, self % numberNeutrons), self % dir(3, self % numberNeutrons))

valeriaRaffuzzi · 2026-02-12T16:16:06Z

ParticleObjects/Source/fileSource_class.f90

+    ! Read and store neutron data from source file
+    ! Value for BroodID is ignored
+    print *, self % numberNeutrons, ' neutrons read from external source file'
+    do i = 1, self % numberNeutrons


Is there a way to avoid looping through the file twice? Does this take long?

I wouldn't think so, to answer both questions. A scan is probably necessary for the allocation.

valeriaRaffuzzi · 2026-02-12T16:18:44Z

ParticleObjects/Source/externalSource_class.f90

+      ! Get pointer to neutron material, return fatal error if not found
+      mat => neutronMaterial_CptrCast(nucData % getMaterial(matIdx))
+      if (.not.associated(mat)) call fatalError(Here, "Nuclear data did not return neutron material.")


The material doesn't seem to be used, don't think we need this.

valeriaRaffuzzi · 2026-02-12T16:19:07Z

ParticleObjects/Source/externalSource_class.f90

+      p % dir = self % dir(:, idx)
+      p % type = P_NEUTRON
+      p % time = ZERO
+      p % wgt = ONE


It would be nice to read the weigth of the particle too!

valeriaRaffuzzi · 2026-02-12T16:20:27Z

ParticleObjects/Source/source_inter.f90

        call dungeon % replace(self % sampleParticle(pRand), i)
      end do
      !$omp end parallel do
+      self % internalCounter = 1


Not sure I understand what this is for, since it never gets updated or incremented?

* Some style changes and code cleaning

valeriaRaffuzzi · 2026-02-13T14:53:46Z

SharedModules/genericProcedures.f90

-    character(100), parameter              :: here = 'readASCII_defReal (genericProcedures.f90)'
+    integer(shortInt), intent(in)                :: unit
+    logical(defBool), intent(in)                 :: readBinary
+    integer(shortInt), intent(out), dimension(:) :: output


Sorry to be pedantic, we tend to do
integer(shortInt), dimension(:), intent(out) :: output

valeriaRaffuzzi · 2026-02-13T14:55:05Z

ParticleObjects/particleDungeon_class.f90

    id = 10
    ! Open the file in requested mode
-    if (writeBinary) then
+    if (writeBinary == 2) then


I like this solution! I think it might be more readable to include a parameter either to this file or to universalVariables like
BINARY_FILE = 2
and then here you could have
if (fileType == BINARY_FILE) then .....

* Added a parameter variable in particleDungeon_class.f90 for writting binary file * Modified CMakeList.txt to hopefully fix the check failure

compilation

ChasingNeutrons · 2026-02-17T13:36:08Z

ParticleObjects/Source/fileSource_class.f90

+  !!   - failed to read source file
+  !!
+  subroutine init(self, dict, geom)
+    class(fileSource), intent(inout)     :: self


Align with the other variables

ChasingNeutrons · 2026-02-17T13:37:39Z

ParticleObjects/Source/fileSource_class.f90

+
+  type, public, extends(source) :: fileSource
+    private
+        integer(shortInt)                            :: numberNeutrons


Probably better as a longInt if we want quite big sources.

ChasingNeutrons · 2026-02-17T13:38:55Z

ParticleObjects/Source/fileSource_class.f90

+    ! Read source data from file, binary or ASCII
+    self % numberNeutrons = 0
+    id = 10
+    ! ASCII file reading


I would add a few more spaces for readability I think

ChasingNeutrons · 2026-02-17T13:44:34Z

ParticleObjects/Source/fileSource_class.f90

+    else
+      nucData => ndReg_getNeutronCE()
+    end if
+    if(.not.associated(nucData)) call fatalError(Here, 'Failed to retrieve Nuclear Database')


It isn't clear to me why knowledge of the nuclear database is necessary for this source. I think it can be removed entirely.

ChasingNeutrons · 2026-02-17T13:47:51Z

ParticleObjects/particleDungeon_class.f90

+
+    id = 10
+    ! Open the file in requested mode
+    if (writeBinary == BINARY_FILE) then


It isn't clear to me why this can't just be a bool, which would remove the need for the parameter defined above.

ChasingNeutrons · 2026-02-17T13:50:57Z

PhysicsPackages/eigenPhysicsPackage_class.f90

+      if (self % printSource /= 0) then
        call self % nextCycle % printToFile(trim(self % outputFile) // '_source' // numToChar(i) // &
-                                            '_rank' // numToChar(getMPIRank()))
+                                            '_rank' // numToChar(getMPIRank()), self % printSource)


Ah, I guess this answers my last question. Maybe replace this with self % printSource == BINARY_FILE, where BINARY_FILE is defined in universalVariables. I think there are a few other places we might want to use that variable in future.

ChasingNeutrons · 2026-02-17T13:51:47Z

PhysicsPackages/eigenPhysicsPackage_class.f90

-    ! Read whether to print particle source per cycle
+    ! Read whether to print particle source per cycle, 1 for ASCII, 2 for binary
    call dict % getOrDefault(self % printSource, 'printSource', 0)
+    if (self % printSource < 0 .or. self % printSource > 2) then


This would likewise use the above mentioned variable.

ChasingNeutrons · 2026-02-17T13:52:13Z

PhysicsPackages/fixedSourcePhysicsPackage_class.f90

-        call self % thisCycle % printToFile(trim(self % outputFile)//'_source'//numToChar(i))
+      ! Print source in ASCII or binary format if requested
+      if (self % printSource /= 0) then
+        call self % thisCycle % printToFile(trim(self % outputFile)//'_source'//numToChar(i), self % printSource)


As for eigen

* Various style fixes * Removed retrieval of dataBase from fileSource

valeriaRaffuzzi

This looks good to me, and if @ChasingNeutrons is also happy with the changes he requested we can merge!

Thanks a lot @bonnettheophile for the development 💯

ChasingNeutrons · 2026-02-17T14:38:32Z

ParticleObjects/Source/fileSource_class.f90

+  use neutronMaterial_inter,   only : neutronMaterial, neutronMaterial_CptrCast
+  use nuclearDataReg_mod,      only : ndReg_getNeutronCE => getNeutronCE, &
+                                      ndReg_getNeutronMG => getNeutronMG
+  use nuclearDatabase_inter,   only : nuclearDatabase


You can also remove references to nuclear data here.

ChasingNeutrons · 2026-02-17T14:39:45Z

ParticleObjects/Source/fileSource_class.f90

+    character(100),parameter :: Here = 'sampleParticle (fileSource_class.f90)'
+
+    ! Sample index of source neutron to be used
+    idx = int(rand % get() * real(self % numberNeutrons, defReal)) + 1


idx needs to be a longInt if there is a large number of particles.

valeriaRaffuzzi requested changes Feb 12, 2026

View reviewed changes

* Renamed externalSource to fileSource

b26bf16

* Some style changes and code cleaning

bonnettheophile requested a review from valeriaRaffuzzi February 13, 2026 14:15

valeriaRaffuzzi requested changes Feb 13, 2026

View reviewed changes

tbonnet and others added 2 commits February 16, 2026 13:57

* Style fix in genericProcedures.f90

da62ce8

* Added a parameter variable in particleDungeon_class.f90 for writting binary file * Modified CMakeList.txt to hopefully fix the check failure

* Modified CMakeList.txt to prevent potential data race upon parallel

87a9e5a

compilation

ChasingNeutrons reviewed Feb 17, 2026

View reviewed changes

* Added parameter variables for source writing in universal variables

9550b94

* Various style fixes * Removed retrieval of dataBase from fileSource

valeriaRaffuzzi approved these changes Feb 17, 2026

View reviewed changes

ChasingNeutrons reviewed Feb 17, 2026

View reviewed changes

valeriaRaffuzzi added 2 commits February 17, 2026 14:45

Remove unnecessary import files in fileSource_class.f90

eff7c63

Change shortInt to longInt for particle number index

1904aad

ChasingNeutrons approved these changes Feb 17, 2026

View reviewed changes

ChasingNeutrons merged commit e8be59d into CambridgeNuclear:main Feb 17, 2026
6 checks passed

		allocate(self % r(3, self % numberNeutrons))
		allocate(self % dir(3, self % numberNeutrons))

Conversation

bonnettheophile commented Feb 9, 2026

Uh oh!

ChasingNeutrons commented Feb 9, 2026

Uh oh!

valeriaRaffuzzi left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

valeriaRaffuzzi left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants