more patches to get the tests running again:

- fixed increment counting in postResults to current output scheme

- corrected node coordinates calculation

- some work on restart (restart test still not running)

code/DAMASK_spectral_driver.f90

@@ -340,6 +340,8 @@ program DAMASK_spectral_Driver
  
 !--------------------------------------------------------------------------------------------------
 ! write header of output file
+ allocate(outputSize(worldsize), source = 0_pInt); outputSize(worldrank+1) = size(materialpoint_results)*8
+ call MPI_Allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
  if (.not. appendToOutFile) then                                                                      ! after restart, append to existing results file
    if (worldrank == 0) then
      open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
@@ -359,20 +361,28 @@ program DAMASK_spectral_Driver
      write(resUnit) 'eoh'    
      close(resUnit)                                                                                   ! end of header
    endif
+   call MPI_File_open(PETSC_COMM_WORLD, &
+                      trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', &
+                      MPI_MODE_WRONLY + MPI_MODE_APPEND, &
+                      MPI_INFO_NULL, &
+                      resUnit, &
+                      ierr)
+   call MPI_File_get_position(resUnit,my_offset,ierr)
+   my_offset = my_offset + sum(outputSize(1:worldrank))
+   call MPI_File_seek (resUnit,my_offset,MPI_SEEK_SET,ierr)
+   call MPI_File_write(resUnit, materialpoint_results, size(materialpoint_results), &
+                       MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
+ else
+   call MPI_File_open(PETSC_COMM_WORLD, &
+                      trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', &
+                      MPI_MODE_WRONLY + MPI_MODE_APPEND, &
+                      MPI_INFO_NULL, &
+                      resUnit, &
+                      ierr)
+   call MPI_File_get_position(resUnit,my_offset,ierr)
+   my_offset = my_offset + sum(outputSize(1:worldrank))
+   call MPI_File_seek (resUnit,my_offset,MPI_SEEK_SET,ierr)
  endif  
- allocate(outputSize(worldsize), source = 0_pInt); outputSize(worldrank+1) = size(materialpoint_results)*8
- call MPI_Allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
- call MPI_File_open(PETSC_COMM_WORLD, &
-                    trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', &
-                    MPI_MODE_WRONLY + MPI_MODE_APPEND, &
-                    MPI_INFO_NULL, &
-                    resUnit, &
-                    ierr)
- call MPI_File_get_position(resUnit,my_offset,ierr)
- my_offset = my_offset + sum(outputSize(1:worldrank))
- call MPI_File_seek (resUnit,my_offset,MPI_SEEK_SET,ierr)
- call MPI_File_write(resUnit, materialpoint_results, size(materialpoint_results), &
-                     MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
  if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0 .and. worldrank == 0_pInt) &
    write(6,'(/,a)') ' header of result file written out'
  flush(6)

code/DAMASK_spectral_solverAL.f90

@@ -218,17 +218,15 @@ subroutine AL_init(temperature)
                                         trim(getSolverJobName()),size(F_lambda_lastInc))
    read (777,rec=1) F_lambda_lastInc
    close (777)
-   if (worldrank == 0_pInt) then
-     call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim))
-     read (777,rec=1) F_aim
-     close (777)
-     call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc))
-     read (777,rec=1) F_aim_lastInc
-     close (777)
-     call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
-     read (777,rec=1) f_aimDot
-     close (777)
-   endif
+   call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim))
+   read (777,rec=1) F_aim
+   close (777)
+   call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc))
+   read (777,rec=1) F_aim_lastInc
+   close (777)
+   call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
+   read (777,rec=1) f_aimDot
+   close (777)
  endif
  
  call utilities_updateIPcoords(F)
@@ -238,8 +236,8 @@ subroutine AL_init(temperature)
  nullify(F_lambda)
  call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)                            ! write data back to PETSc
                              
- if (restartInc > 1_pInt .and. worldrank == 0_pInt) then                                            ! using old values from files
-   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
+ if (restartInc > 1_pInt) then                                                                      ! using old values from files
+   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
      write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
      'reading more values of increment', restartInc - 1_pInt, 'from file'
    flush(6)

code/DAMASK_spectral_solverBasicPETSc.f90

@@ -193,11 +193,9 @@ subroutine basicPETSc_init(temperature)
    call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
    read (777,rec=1) F_lastInc
    close (777)
-   if (worldrank == 0_pInt) then
-     call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
-     read (777,rec=1) f_aimDot
-     close (777)
-   endif
+   call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
+   read (777,rec=1) f_aimDot
+   close (777)
    F_aim         = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3])                         ! average of F
    F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt                                 ! average of F_lastInc 
  endif
@@ -213,8 +211,8 @@ subroutine basicPETSc_init(temperature)
     math_I3)
  call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                 ! write data back to PETSc
 
- if (restartInc > 1_pInt .and. worldrank == 0_pInt) then                                                                      ! using old values from files                                                    
-   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
+ if (restartInc > 1_pInt) then                                                                      ! using old values from files                                                    
+   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
      write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
      'reading more values of increment', restartInc - 1_pInt, 'from file'
    flush(6)

code/DAMASK_spectral_solverPolarisation.f90

@@ -217,17 +217,15 @@ subroutine Polarisation_init(temperature)
                                         trim(getSolverJobName()),size(F_tau_lastInc))
    read (777,rec=1) F_tau_lastInc
    close (777)
-   if (worldrank == 0_pInt) then
-     call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim))
-     read (777,rec=1) F_aim
-     close (777)
-     call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc))
-     read (777,rec=1) F_aim_lastInc
-     close (777)
-     call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
-     read (777,rec=1) f_aimDot
-     close (777)
-   endif
+   call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim))
+   read (777,rec=1) F_aim
+   close (777)
+   call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc))
+   read (777,rec=1) F_aim_lastInc
+   close (777)
+   call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
+   read (777,rec=1) f_aimDot
+   close (777)
  endif
 
  call Utilities_updateIPcoords(F)
@@ -237,8 +235,8 @@ subroutine Polarisation_init(temperature)
  nullify(F_tau)
  call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)                            ! write data back to PETSc
 
- if (restartInc > 1_pInt .and. worldrank == 0_pInt) then                                            ! using old values from files
-   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
+ if (restartInc > 1_pInt) then                                            ! using old values from files
+   if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
      write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
      'reading more values of increment', restartInc - 1_pInt, 'from file'
    flush(6)

code/mesh.f90

@@ -1315,18 +1315,23 @@ subroutine mesh_spectral_build_nodes(fileUnit)
 
  implicit none
  integer(pInt),              intent(in) :: fileUnit
- integer(pInt)                          :: n, i, j, k
+ integer(pInt)                          :: n
 
  allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal)
  allocate (mesh_node  (3,mesh_Nnodes), source = 0.0_pReal)
  
- n = 0_pInt
- do k = 1, gridLocal(3); do j = 1, gridLocal(2); do i = 1, gridLocal(1)
-   n = n + 1_pInt
-   mesh_node0(1,n) = real(i)*geomSizeLocal(1)/real(gridLocal(1))
-   mesh_node0(2,n) = real(j)*geomSizeLocal(2)/real(gridLocal(2))
-   mesh_node0(3,n) = real(k)*geomSizeLocal(3)/real(gridLocal(3)) + geomSizeOffset
- enddo; enddo; enddo 
+ forall (n = 0_pInt:mesh_Nnodes-1_pInt)
+   mesh_node0(1,n+1_pInt) = mesh_unitlength * &
+           geomSizeLocal(1)*real(mod(n,(gridLocal(1)+1_pInt) ),pReal) &
+                                                  / real(gridLocal(1),pReal)
+   mesh_node0(2,n+1_pInt) = mesh_unitlength * &
+           geomSizeLocal(2)*real(mod(n/(gridLocal(1)+1_pInt),(gridLocal(2)+1_pInt)),pReal) &
+                                                  / real(gridLocal(2),pReal)
+   mesh_node0(3,n+1_pInt) = mesh_unitlength * &
+           geomSizeLocal(3)*real(mod(n/(gridLocal(1)+1_pInt)/(gridLocal(2)+1_pInt),(gridLocal(3)+1_pInt)),pReal) &
+                                                  / real(gridLocal(3),pReal) + &
+           geomSizeOffset                                       
+ end forall 
 
  mesh_node = mesh_node0
 
@@ -1599,14 +1604,7 @@ function mesh_regrid(adaptive,resNewInput,minRes)
 !--------------------------------------------------------------------------------------------------
 ! read in deformation gradient to calculate coordinates, shape depend of selected solver
  select case(spectral_solver)
-   case('basic')
-     allocate(spectralF33(3,3,grid(1),grid(2),grid(3)))
-     call IO_read_realFile(FILEUNIT,'F',trim(getSolverJobName()),size(spectralF33))
-     read (FILEUNIT,rec=1) spectralF33
-     close (FILEUNIT)
-     Favg = sum(sum(sum(spectralF33,dim=5),dim=4),dim=3) * wgt
-     coordinates = reshape(mesh_deformedCoordsFFT(geomSize,spectralF33),[3,mesh_NcpElems])
-   case('basicpetsc','al')
+   case('basicpetsc','al','polarization')
      allocate(spectralF9(9,grid(1),grid(2),grid(3)))
      call IO_read_realFile(FILEUNIT,'F',trim(getSolverJobName()),size(spectralF9))
      read (FILEUNIT,rec=1) spectralF9
@@ -1750,15 +1748,8 @@ function mesh_regrid(adaptive,resNewInput,minRes)
  
 !--------------------------------------------------------------------------------------------------
 ! set F to average values
- select case(spectral_solver)
-   case('basic')
-     allocate(spectralF33New(3,3,resNew(1),resNew(2),resNew(3)))
-     spectralF33New = spread(spread(spread(Favg,3,resNew(1)),4,resNew(2)),5,resNew(3))
-     call IO_write_jobRealFile(FILEUNIT,'F',size(spectralF33New))
-     write (FILEUNIT,rec=1) spectralF33New
-     close (FILEUNIT)
-     
-   case('basicpetsc','al')
+ select case(spectral_solver)     
+   case('basicpetsc','al','polarization')
      allocate(spectralF9New(9,resNew(1),resNew(2),resNew(3)))
      spectralF9New = spread(spread(spread(reshape(Favg,[9]),2,resNew(1)),3,resNew(2)),4,resNew(3))
      call IO_write_jobRealFile(FILEUNIT,'F',size(spectralF9New))

processing/post/postResults.py

@@ -138,7 +138,7 @@ class MPIEspectral_result:    # mimic py_post result object
     if self.element_scalars == None:
       self.N_element_scalars = self._keyedPackedArray('materialpoint_sizeResults',count=1,type='i',default=0)[0]
 
-    self.N_positions  = (self.filesize-self.dataOffset)/(8+self.N_elements*self.N_element_scalars*8)
+    self.N_positions  = (self.filesize-self.dataOffset)/(self.N_elements*self.N_element_scalars*8)
     self.N_increments = 1                                                    # add zero'th entry
     for i in range(self.N_loadcases):
       self.N_increments += self._increments[i]//self._frequencies[i]
@@ -268,7 +268,7 @@ class MPIEspectral_result:    # mimic py_post result object
 
   def element_scalar(self,e,idx):
     incStart =  self.dataOffset \
-             +  self.position*8*(self.N_elements*self.N_element_scalars)
+             +  self.position*8*self.N_elements*self.N_element_scalars
                                 # header & footer + extra header and footer for 4 byte int range (Fortran)
                                 # values
     where = (e*self.N_element_scalars + idx)*8