corrected determination of restartRead from commandline (now done by FEsolving instead of DAMASK_spectral_interface)

recompiled fftw libraries with --enable-shared to make them work with f2py
added draft for compilation test of spectral solver

code/DAMASK_spectral.f90

@@ -68,9 +68,10 @@ program DAMASK_spectral
  integer(pInt), dimension (1 + maxNchunksLoadcase*2) :: posLoadcase
  integer(pInt), parameter :: maxNchunksGeom = 7_pInt                                  ! 4 identifiers, 3 values
  integer(pInt), dimension (1 + maxNchunksGeom*2) ::  posGeom
- integer(pInt) :: myUnit, N_l, N_s, N_t, N_n, N_Fdot, headerLength                    ! numbers of identifiers
+ integer(pInt) :: headerLength,N_l=0_pInt, N_t=0_pInt, N_n=0_pInt, N_Fdot=0_pInt
+ integer(pInt), parameter :: myUnit = 234_pInt
  character(len=1024) :: path, line, keyword
- logical ::  gotResolution, gotDimension, gotHomogenization
+ logical ::  gotResolution =.false., gotDimension =.false., gotHomogenization = .false.
  
 ! variables storing information from loadcase file
 !ToDo: create Data Structure loadcase
@@ -91,14 +92,14 @@ program DAMASK_spectral
 
 ! variables storing information from geom file
  real(pReal) :: wgt
- real(pReal), dimension(3) ::  geomdimension    ! physical dimension of volume element in each direction
+ real(pReal), dimension(3) ::  geomdimension = 0.0_pReal                    ! physical dimension of volume element in each direction
  integer(pInt) :: homog                                                     ! homogenization scheme used
- integer(pInt), dimension(3) :: resolution      ! resolution (number of Fourier points) in each direction
+ integer(pInt), dimension(3) :: resolution = 1_pInt                         ! resolution (number of Fourier points) in each direction
- logical :: spectralPictureMode                 ! indicating 1 to 1 mapping of FP to microstructure
+ logical :: spectralPictureMode = .false.                                   ! indicating 1 to 1 mapping of FP to microstructure
 
-! stress etc.
+! stress, stiffness and compliance average etc.
  real(pReal), dimension(3,3) ::                pstress, pstress_av, defgrad_av, &
-                                                        defgradAim, defgradAimOld, defgradAimCorr,&
+                                               defgradAim = math_I3, defgradAimOld= math_I3, defgradAimCorr= math_I3,&
                                                mask_stress, mask_defgrad, fDot, &
                                                pstress_av_load, defgradAim_lab               ! quantities rotated to other coordinate system
  real(pReal), dimension(3,3,3,3) ::            dPdF, c0_reference, c_current, s_prev, c_prev ! stiffness and compliance
@@ -106,7 +107,7 @@ program DAMASK_spectral
  real(pReal), dimension(6,6) ::                dsde                                          ! small strain stiffness
  real(pReal), dimension(9,9) ::                s_prev99, c_prev99                            ! compliance and stiffness in matrix notation 
  real(pReal), dimension(:,:), allocatable ::   s_reduced, c_reduced                          ! reduced compliance and stiffness (only for stress BC)
- integer(pInt) ::                                       size_reduced                                  ! number of stress BCs
+ integer(pInt) ::                              size_reduced = 0.0_pReal                                  ! number of stress BCs
 
 ! pointwise data 
  real(pReal), dimension(:,:,:,:,:), allocatable ::      workfft, defgrad, defgradold
@@ -150,15 +151,8 @@ program DAMASK_spectral
  !$OMP END CRITICAL (write2out)
 
 ! Reading the loadcase file and allocate variables
- myUnit = 234_pInt
  path = getLoadcaseName()
  if (.not. IO_open_file(myUnit,path)) call IO_error(error_ID=30,ext_msg = trim(path))
- 
- N_l = 0_pInt
- N_Fdot = 0_pInt
- N_t = 0_pInt
- N_n = 0_pInt
-
  rewind(myUnit)
  do
    read(myUnit,'(a1024)',END = 100) line
@@ -261,17 +255,10 @@ program DAMASK_spectral
 101 close(myUnit)
 
 !read header of geom file to get the information needed before the complete geom file is intepretated by mesh.f90
- gotResolution =.false.
- gotDimension =.false.
- gotHomogenization = .false.
- spectralPictureMode = .false.
- resolution = 1_pInt
- geomdimension = 0.0_pReal
  
  path = getModelName()
  if (.not. IO_open_file(myUnit,trim(path)//InputFileExtension))&
         call IO_error(error_ID=101,ext_msg = trim(path)//InputFileExtension)
-
  rewind(myUnit)
  read(myUnit,'(a1024)') line
  posGeom = IO_stringPos(line,2)
@@ -332,9 +319,6 @@ program DAMASK_spectral
  allocate (xi         (3,resolution(1)/2+1,resolution(2),resolution(3)));  xi          =0.0_pReal 
 
  wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal)
- defgradAim    = math_I3
- defgradAimOld = math_I3
- defgrad_av    = math_I3
  
 ! Initialization of CPFEM_general (= constitutive law) and of deformation gradient field
  call CPFEM_initAll(bc_temperature(1),1_pInt,1_pInt)

code/DAMASK_spectral_interface.f90

@@ -27,9 +27,7 @@ MODULE DAMASK_interface
  character(len=64), parameter :: FEsolver = 'Spectral'
  character(len=5),  parameter :: InputFileExtension = '.geom'
  character(len=4),  parameter :: LogFileExtension = '.log'    !until now, we don't have a log file. But IO.f90 requires it
- logical :: restartReadFlag
  character(len=1024) :: geometryParameter,loadcaseParameter
- integer(pInt) :: restartReadStep
 CONTAINS
 
 !********************************************************************
@@ -41,11 +39,7 @@ subroutine DAMASK_interface_init()
  implicit none
 
  character(len=1024) commandLine
- integer(pInt):: i, start, length
+ integer(pInt):: i, start = 0_pInt, length=0_pInt
-
- start = 0_pInt
- length= 0_pInt
- restartReadFlag = .false.
 
  call get_command(commandLine)
 
@@ -66,7 +60,6 @@ subroutine DAMASK_interface_init()
  geometryParameter = ''                                            ! should be empty
  geometryParameter(1:length)=commandLine(start:start+length)
  
- call get_command(commandLine)
  do i=1,len(commandLine)                                           ! remove capitals
    if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32)
  enddo
@@ -85,24 +78,6 @@ subroutine DAMASK_interface_init()
  loadcaseParameter = ''                                            ! should be empty
  loadcaseParameter(1:length)=commandLine(start:start+length)
  
- do i=1,len(commandLine)                                           ! remove capitals
-   if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32)
- enddo
- 
- start = index(commandLine,'-r',.true.) + 3_pInt                   ! search for '-r' and jump forward to given name
- if (index(commandLine,'--restart',.true.)>0) then                 ! if '--restart' is found, use that (contains '-r')
-   start = index(commandLine,'--restart',.true.) + 10_pInt
- endif               
- length = index(commandLine(start:len(commandLine)),' ',.false.)
-
- restartReadStep = 0_pInt
- if(start/=3_pInt) then                                           ! found -r
-   read(commandLine(start:start+length),'(I12)') restartReadStep
-   restartReadFlag = .true.
- endif
- 
- if(restartReadStep<1_pInt .and. RestartReadFlag .eq. .true.) stop 'Invalid Restart Parameter, terminating DAMASK' ! Functions from IO.f90 are not available
-
  !$OMP CRITICAL (write2out)
  write(6,*)
  write(6,*) '<<<+-  DAMASK_spectral_interface init  -+>>>'
@@ -110,11 +85,6 @@ subroutine DAMASK_interface_init()
  write(6,*)
  write(6,*) 'Geometry Parameter: ', trim(geometryParameter)
  write(6,*) 'Loadcase Parameter: ', trim(loadcaseParameter)
- if (restartReadFlag) then
-   write(6,*) 'Restart Read: ', restartReadFlag
- else 
-   write(6,'(a,I5)') ' Restart Read at Step: ', restartReadStep
- endif
  write(6,*)
  !$OMP END CRITICAL (write2out)
 

code/FEsolving.f90

@@ -25,7 +25,7 @@
  use prec, only: pInt,pReal
  implicit none
 
- integer(pInt) :: cycleCounter = 0_pInt, theInc = -1_pInt
+ integer(pInt) :: cycleCounter = 0_pInt, theInc = -1_pInt, restartReadStep = 0_pInt
  real(pReal)   :: theTime = 0.0_pReal, theDelta = 0.0_pReal
  logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.,terminallyIll = .false.
  logical :: symmetricSolver = .false. 
@@ -54,16 +54,29 @@
  
  integer(pInt), parameter :: fileunit = 222
  integer(pInt), parameter :: maxNchunks = 6
- integer(pInt) i
+ integer(pInt):: i, start = 0_pInt, length=0_pInt
  integer(pInt), dimension(1+2*maxNchunks) :: positions
  character(len=64) tag
- character(len=1024) line
+ character(len=1024) line, commandLine
 
  FEmodelGeometry = getModelName()
  if (IO_open_inputFile(fileunit,FEmodelGeometry)) then
    if(trim(FEsolver)=='Spectral') then
      restartWrite = .true.
-     restartRead  = restartReadFlag
+     call get_command(commandLine)                                     ! may contain capitals
+     do i=1,len(commandLine)                                           ! remove capitals
+       if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32)
+     enddo
+     start = index(commandLine,'-r',.true.) + 3_pInt                   ! search for '-r' and jump forward to given name
+     if (index(commandLine,'--restart',.true.)>0) then                 ! if '--restart' is found, use that (contains '-r')
+       start = index(commandLine,'--restart',.true.) + 10_pInt
+     endif               
+     length = index(commandLine(start:len(commandLine)),' ',.false.)
+     if(start/=3_pInt) then                                            ! found at least -r
+       read(commandLine(start:start+length),'(I12)') restartReadStep
+       restartRead = .true.
+     endif
+   if(restartReadStep<0_pInt .and. restartRead .eq. .true.)  call IO_error(error_ID=47)
    else
      rewind(fileunit)
      do

code/IO.f90

@@ -1166,6 +1166,8 @@ endfunction
    msg = 'incomplete information in spectral mesh header'
  case (46)
    msg = 'not a rotation defined for loadcase rotation'
+ case (47)
+   msg = 'invalid restart increment given'
  case (50)
    msg = 'writing constitutive output description'
  case (100)

code/makefile

@@ -6,8 +6,8 @@
 # Uses OpenMP to parallelise the material subroutines (set number of cores with "export DAMASK_NUM_THREADS=n" to n)
 # Uses linux threads to parallelise fftw3
 # Besides of the f90 files written at MPIE, the two library files of fftw3 "libfftw3_threads.a" "libfftw3.a" are also needed
-# Install fftw3 (v3.3 is tested) with "./configure --enable-threads --enable-sse2" and "make"; "make install" is not needed
+# Install fftw3 (v3.3 is tested) with "./configure --enable-threads --enable-sse2 --enable-shared" and "make"; "make install" is not needed
-# as long as the two library files "libfftw3_threads.a" "libfftw3.a" are copied to the code/include directory.
+# as long as the two library files "libfftw3_threads.a" "libfftw3.a" are copied to the /../lib directory.
 # Need the AMD Core Math Library to be installed (v 4.4 is tested)
 ########################################################################################
 # OPTIONS = standard (alternative): meaning

processing/post/make_postprocessingMath

@@ -10,8 +10,7 @@
 # for the generation of the pyf file:
 # f2py -m postprocessingMath -h postprocessingMath.pyf --overwrite-signature postprocessingMath.f90
 
-# --f90flags="-heap-arrays 500000000" prevents segmentation fault for large arrays
+# use ./configure --enable-sse2 --enable-shared for the compilation of fftw3.3
-# use ./configure --enable-portable-binary --enable-shared for the compilation of fftw3.2.2
 
 if [[ $# -eq 0 ]]; then
   wd='.'
@@ -20,10 +19,10 @@ else
 fi
 
 cd $wd
-rm postprocessingMath.so
+rm ../../lib/postprocessingMath.so
 f2py -c \
 postprocessingMath.pyf \
 postprocessingMath.f90 \
--lfftw3 \
+../../lib/libfftw3.a \
 -L./
-#--f90flags="-heap-arrays 500000000"
+mv postprocessingMath.so ../../lib/.