restructured algorithm: moved into loop to reallocate fields and replan FFTW in case resolution changes during runtime ==> regridding

introduced parameters for selective debugging of spectral code and partly introduced the advanced divergence calculation again which is controlled by debug.config
added switch in numerics to control divergence behavior (uncorrected and corrected by phenomenological factor)
added precision directive to all values I found

code/IO.f90

@@ -1416,6 +1416,8 @@ endfunction
  select case (warning_ID)
  case (33_pInt)
    msg = 'cannot guess along trajectory for first step of first loadcase'
+ case (47_pInt)
+   msg = 'No valid parameter for FFTW given, using FFTW_PATIENT'
  case (101_pInt)
    msg = '+    crystallite debugging off...    +'
  case (600_pInt)

code/config/debug.config

@@ -15,3 +15,7 @@ selective        1                    # >0 true to switch on e,i,g selective deb
 element          1                    # selected element for debugging (synonymous: "el", "e")
 ip               1                    # selected integration point for debugging (synonymous: "integrationpoint", "i")
 grain            1                    # selected grain at ip for debugging (synonymous: "gr", "g")
+
+### spectral solver debugging parameters ###
+generalDebugSpectral    0             # > 0: general (algorithmical) debug outputs
+divergenceDebugSpectral 0             # > 0: calculate more divergence measures and print them out

code/debug.f90

@@ -51,6 +51,7 @@ real(pReal) :: debug_jacobianMax
 real(pReal) :: debug_jacobianMin
 logical :: debug_selectiveDebugger = .true.
 integer(pInt) :: debug_verbosity = 1_pInt
+integer(pInt) :: spectral_debug_verbosity = 0_pInt
 
 CONTAINS
 
@@ -122,6 +123,12 @@ subroutine debug_init()
               debug_selectiveDebugger = IO_intValue(line,positions,2) > 0_pInt
         case ('verbosity')
               debug_verbosity = IO_intValue(line,positions,2)
+        case ('generaldebugspectral')                          ! use bitwise logical and, continue with +8_pInt
+              if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 1_pInt 
+        case ('divergencedebugspectral')
+              if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 2_pInt
+        case ('restartdebugspectral')
+              if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 4_pInt
       endselect
     enddo
     100 close(fileunit)
@@ -164,7 +171,11 @@ subroutine debug_init()
     debug_i = 0_pInt
     debug_g = 0_pInt
   endif
-
+  !$OMP CRITICAL (write2out)                                                   ! bitwise coded
+  if (iand(spectral_debug_verbosity,1_pInt)==1_pInt)  write(6,'(a)') ' Spectral General Debugging'
+  if (iand(spectral_debug_verbosity,2_pInt)==2_pInt)  write(6,'(a)') ' Spectral Divergence Debugging'
+  if (iand(spectral_debug_verbosity,4_pInt)==4_pInt)  write(6,'(a)') ' Spectral Restart Debugging'
+  !$OMP END CRITICAL (write2out)
 
 endsubroutine
  
@@ -196,6 +207,7 @@ subroutine debug_reset()
   debug_stressMin = huge(1.0_pReal)
   debug_jacobianMax = -huge(1.0_pReal)
   debug_jacobianMin = huge(1.0_pReal)
+  spectral_debug_verbosity = 0.0_pReal
 
 
 endsubroutine

code/numerics.f90

@@ -70,9 +70,10 @@ real(pReal) ::                  relevantStrain, &                       ! strain
                                 err_stress_tolrel, &                    ! factor to multiply with highest stress to get err_stress_tol
                                 fftw_timelimit, &                       ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org
                                 rotation_tol                            ! tolerance of rotation specified in loadcase
-character(len=64)               fftw_planner_flag                       ! sets the planig-rigor flag, see manual on www.fftw.org
+character(len=64) ::            fftw_planner_flag                       ! sets the planig-rigor flag, see manual on www.fftw.org
-logical                         memory_efficient                        ! for fast execution (pre calculation of gamma_hat)
+logical ::                      memory_efficient,&                      ! for fast execution (pre calculation of gamma_hat)
-integer(pInt)                   itmax , &                               ! maximum number of iterations
+                                divergence_correction                     ! correct divergence calculation in fourier space
+integer(pInt) ::                itmax , &                               ! maximum number of iterations
 
 
 !* Random seeding parameters
@@ -172,8 +173,8 @@ subroutine numerics_init()
   fftw_timelimit          = -1.0_pReal   ! no timelimit of plan creation for FFTW
   fftw_planner_flag       ='FFTW_PATIENT'
   rotation_tol            = 1.0e-12 
-
+  divergence_correction   = .true.
-!* Random seeding parameters: added <<<updated 27.08.2009>>>
+!* Random seeding parameters
   fixedSeed               = 0_pInt
 
 
@@ -290,6 +291,8 @@ subroutine numerics_init()
               fftw_planner_flag = IO_stringValue(line,positions,2)
         case ('rotation_tol')
               rotation_tol = IO_floatValue(line,positions,2)
+        case ('divergence_correction')
+              divergence_correction = IO_intValue(line,positions,2)  > 0_pInt
 
 !* Random seeding parameters
         case ('fixed_seed')
@@ -326,15 +329,13 @@ subroutine numerics_init()
     write(6,'(a24,x,e8.1)') ' rTol_crystalliteTemp:   ',rTol_crystalliteTemperature
     write(6,'(a24,x,e8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
     write(6,'(a24,x,e8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
-    write(6,'(a24,2(x,i8))')'integrator:             ',numerics_integrator
+    write(6,'(a24,2(x,i8),/)')' integrator:             ',numerics_integrator
-    write(6,*)
   
     write(6,'(a24,x,i8)')   ' nHomog:                 ',nHomog
     write(6,'(a24,x,e8.1)') ' subStepMinHomog:        ',subStepMinHomog
     write(6,'(a24,x,e8.1)') ' subStepSizeHomog:       ',subStepSizeHomog
     write(6,'(a24,x,e8.1)') ' stepIncreaseHomog:      ',stepIncreaseHomog
-    write(6,'(a24,x,i8)')   'nMPstate:               ',nMPstate
+    write(6,'(a24,x,i8,/)') ' nMPstate:               ',nMPstate
-    write(6,*)
 
 !* RGC parameters
     write(6,'(a24,x,e8.1)') ' aTol_RGC:               ',absTol_RGC
@@ -348,8 +349,7 @@ subroutine numerics_init()
     write(6,'(a24,x,e8.1)') ' maxrelaxation_RGC:      ',maxdRelax_RGC
     write(6,'(a24,x,e8.1)') ' maxVolDiscrepancy_RGC:  ',maxVolDiscr_RGC
     write(6,'(a24,x,e8.1)') ' volDiscrepancyMod_RGC:  ',volDiscrMod_RGC
-    write(6,'(a24,x,e8.1)') 'discrepancyPower_RGC:   ',volDiscrPow_RGC
+    write(6,'(a24,x,e8.1,/)') ' discrepancyPower_RGC:   ',volDiscrPow_RGC
-    write(6,*)
 
 !* spectral parameters
     write(6,'(a24,x,e8.1)')   ' err_div_tol:            ',err_div_tol
@@ -363,16 +363,14 @@ subroutine numerics_init()
     endif
     write(6,'(a24,x,a)')      ' fftw_planner_flag:      ',trim(fftw_planner_flag)
     write(6,'(a24,x,e8.1)')   ' rotation_tol:           ',rotation_tol
-    write(6,*)
+    write(6,'(a24,x,L8,/)')   ' divergence_correction:  ',divergence_correction
 
 !* Random seeding parameters
-    write(6,'(a24,x,i16)')   'fixed_seed:             ',fixedSeed
+    write(6,'(a24,x,i16,/)')   ' fixed_seed:             ',fixedSeed
-    write(6,*)
   !$OMP END CRITICAL (write2out)
 
 !* openMP parameter
-!$  write(6,'(a24,x,i8)')   'number of threads:      ',DAMASK_NumThreadsInt
+!$  write(6,'(a24,x,i8,/)')   ' number of threads:      ',DAMASK_NumThreadsInt
-!$  write(6,*)
   
   ! sanity check  
   if (relevantStrain <= 0.0_pReal)          call IO_error(260)