From e01d271ee4c668fc87c58cfcebf3e8f3d35676c1 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 27 Jul 2021 07:43:35 +0200
Subject: [PATCH] store geometry (for full reproducibility)

---
 src/grid/discretization_grid.f90 | 39 ++++++++++++++++++++------------
 1 file changed, 24 insertions(+), 15 deletions(-)

diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90
index 049b3b60d..f94c9921b 100644
--- a/src/grid/discretization_grid.f90
+++ b/src/grid/discretization_grid.f90
@@ -68,11 +68,28 @@ subroutine discretization_grid_init(restart)
     devNull, z, z_offset
   integer, dimension(worldsize) :: &
     displs, sendcounts
+  integer :: fileUnit, myStat
+  integer(pI64) :: &
+    fileLength
+  character(len=:), allocatable :: &
+    fileContent
+
 
   print'(/,a)', ' <<<+-  discretization_grid init  -+>>>'; flush(IO_STDOUT)
 
+
   if(worldrank == 0) then
-    call readVTI(grid,geomSize,origin,materialAt_global)
+    inquire(file = trim(interface_geomFile), size=fileLength)
+    open(newunit=fileUnit, file=trim(interface_geomFile), access='stream',&
+         status='old', position='rewind', action='read',iostat=myStat)
+    if(myStat /= 0) call IO_error(100,ext_msg=trim(interface_geomFile))
+    allocate(character(len=fileLength)::fileContent)
+    read(fileUnit) fileContent
+    close(fileUnit)
+    call readVTI(grid,geomSize,origin,materialAt_global,fileContent)
+    call results_openJobFile(parallel=.false.)
+    call results_writeDataset_str(fileContent,'setup',interface_geomFile,'geometry definition (grid solver)')
+    call results_closeJobFile
   else
     allocate(materialAt_global(0))                                                                  ! needed for IntelMPI
   endif
@@ -157,7 +174,8 @@ end subroutine discretization_grid_init
 !> @brief Parse vtk image data (.vti)
 !> @details https://vtk.org/Wiki/VTK_XML_Formats
 !--------------------------------------------------------------------------------------------------
-subroutine readVTI(grid,geomSize,origin,material)
+subroutine readVTI(grid,geomSize,origin,material, &
+                   fileContent)
 
   integer,     dimension(3), intent(out) :: &
     grid                                                                                            ! grid   (across all processes!)
@@ -166,28 +184,19 @@ subroutine readVTI(grid,geomSize,origin,material)
     origin                                                                                          ! origin (across all processes!)
   integer,     dimension(:), intent(out), allocatable :: &
     material
+  character(len=*),          intent(in) :: &
+    fileContent
 
-  character(len=:), allocatable :: fileContent, dataType, headerType
+  character(len=:), allocatable :: dataType, headerType
   logical :: inFile,inImage,gotCellData,compressed
-  integer :: fileUnit, myStat
   integer(pI64) :: &
-    fileLength, &                                                                                   !< length of the geom file (in characters)
     startPos, endPos, &
     s
 
+
   grid = -1
   geomSize = -1.0_pReal
 
-!--------------------------------------------------------------------------------------------------
-! read raw data as stream
-  inquire(file = trim(interface_geomFile), size=fileLength)
-  open(newunit=fileUnit, file=trim(interface_geomFile), access='stream',&
-       status='old', position='rewind', action='read',iostat=myStat)
-  if(myStat /= 0) call IO_error(100,ext_msg=trim(interface_geomFile))
-  allocate(character(len=fileLength)::fileContent)
-  read(fileUnit) fileContent
-  close(fileUnit)
-
   inFile         = .false.
   inImage        = .false.
   gotCelldata    = .false.