Merge branch '39-simplify-obscure-numerics-integration-mode' into 'development'

Resolve "simplify obscure numerics integration mode"

Closes #39

See merge request damask/DAMASK!36

src/crystallite.f90

@@ -114,6 +114,7 @@ module crystallite
  end enum
  integer(kind(undefined_ID)),dimension(:,:),   allocatable,          private :: &
    crystallite_outputID                                                                             !< ID of each post result output
+ procedure(), pointer :: integrateState
 
  public :: &
    crystallite_init, &
@@ -122,6 +123,7 @@ module crystallite
    crystallite_push33ToRef, &
    crystallite_postResults
  private :: &
+   integrateState, &
    crystallite_integrateStateFPI, &
    crystallite_integrateStateEuler, &
    crystallite_integrateStateAdaptiveEuler, &
@@ -149,6 +151,7 @@ subroutine crystallite_init
    debug_crystallite, &
    debug_levelBasic
  use numerics, only: &
+   numerics_integrator, &
    worldrank, &
    usePingPong
  use math, only: &
@@ -271,6 +274,20 @@ subroutine crystallite_init
  allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), &
                                      size(config_crystallite)), source=0_pInt)
 
+ select case(numerics_integrator(1))
+   case(1_pInt)
+     integrateState => crystallite_integrateStateFPI
+   case(2_pInt)
+     integrateState => crystallite_integrateStateEuler
+   case(3_pInt)
+     integrateState => crystallite_integrateStateAdaptiveEuler
+   case(4_pInt)
+     integrateState => crystallite_integrateStateRK4
+   case(5_pInt)
+     integrateState => crystallite_integrateStateRKCK45
+ end select
+
+
 
  do c = 1_pInt, size(config_crystallite)
 #if defined(__GFORTRAN__)
@@ -509,8 +526,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
    subStepMinCryst, &
    subStepSizeCryst, &
    stepIncreaseCryst, &
-   numerics_integrator, &
-   numerics_integrationMode, &
    numerics_timeSyncing
  use debug, only: &
    debug_level, &
@@ -663,7 +678,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
  endif singleRun
 
  NiterationCrystallite = 0_pInt
- numerics_integrationMode = 1_pInt
  cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2))))
 
    if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
@@ -1041,25 +1055,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
 
    ! --- integrate --- requires fully defined state array (basic + dependent state)
 
-   if (any(crystallite_todo)) then
+   if (any(crystallite_todo)) call integrateState()
-     if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then
-       write(6,'(/,a,i3)') '<< CRYST >> using state integrator ',numerics_integrator(numerics_integrationMode)
-       flush(6)
-     endif
-     select case(numerics_integrator(numerics_integrationMode))
-       case(1_pInt)
-         call crystallite_integrateStateFPI()
-       case(2_pInt)
-         call crystallite_integrateStateEuler()
-       case(3_pInt)
-         call crystallite_integrateStateAdaptiveEuler()
-       case(4_pInt)
-         call crystallite_integrateStateRK4()
-       case(5_pInt)
-         call crystallite_integrateStateRKCK45()
-     end select
-   endif
-
    where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) &                  ! do not try non-converged & fully cutbacked any further
      crystallite_todo = .true.
 
@@ -2025,8 +2021,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
    debug_levelExtensive, &
    debug_levelSelective
  use numerics, only: &
-   rTol_crystalliteState, &
+   rTol_crystalliteState
-   numerics_integrationMode
  use FEsolving, only: &
    FEsolving_execElem, &
    FEsolving_execIP
@@ -2094,7 +2089,6 @@ subroutine crystallite_integrateStateAdaptiveEuler()
  sourceStateResiduum = 0.0_pReal
  relSourceStateResiduum = 0.0_pReal
 
- integrationMode: if (numerics_integrationMode == 1_pInt) then
 
  !$OMP PARALLEL
    ! --- DOT STATE (EULER INTEGRATION) ---
@@ -2194,7 +2188,6 @@ subroutine crystallite_integrateStateAdaptiveEuler()
      enddo; enddo; enddo
    !$OMP ENDDO
  !$OMP END PARALLEL
- endif integrationMode
 
 
  ! --- STRESS INTEGRATION (EULER INTEGRATION) ---
@@ -2214,9 +2207,6 @@ subroutine crystallite_integrateStateAdaptiveEuler()
    enddo; enddo; enddo
  !$OMP END PARALLEL DO
 
-
- if (numerics_integrationMode == 1_pInt) then
-
    !$OMP PARALLEL
    ! --- DOT STATE (HEUN METHOD) ---
 
@@ -2335,17 +2325,6 @@ subroutine crystallite_integrateStateAdaptiveEuler()
    !$OMP ENDDO
  !$OMP END PARALLEL
 
- elseif (numerics_integrationMode > 1) then ! stiffness calculation
-
-   !$OMP PARALLEL DO
-     do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e)                    ! iterate over elements, ips and grains
-       crystallite_converged(g,i,e) = crystallite_todo(g,i,e) .or. crystallite_converged(g,i,e)              ! ... converged per definitionem
-     enddo; enddo; enddo
-   !$OMP END PARALLEL DO
-
- endif
-
-
 
  ! --- NONLOCAL CONVERGENCE CHECK ---
 
@@ -2376,7 +2355,6 @@ subroutine crystallite_integrateStateEuler()
    debug_levelExtensive, &
    debug_levelSelective
  use numerics, only: &
-   numerics_integrationMode, &
    numerics_timeSyncing
  use FEsolving, only: &
    FEsolving_execElem, &
@@ -2423,7 +2401,6 @@ eIter = FEsolving_execElem(1:2)
 
  singleRun = (eIter(1) == eIter(2) .and. iIter(1,eIter(1)) == iIter(2,eIter(2)))
 
- if (numerics_integrationMode == 1_pInt) then
  !$OMP PARALLEL
 
    ! --- DOT STATE  ---
@@ -2528,7 +2505,6 @@ eIter = FEsolving_execElem(1:2)
    enddo; enddo; enddo
    !$OMP ENDDO
   !$OMP END PARALLEL
- endif
 
 
  !$OMP PARALLEL

src/numerics.f90

@@ -26,9 +26,8 @@ module numerics
    worldsize                  =  0_pInt                                                             !< MPI worldsize (/=0 for MPI simulations only)
  integer(4), protected, public :: &
    DAMASK_NumThreadsInt       =  0                                                                  !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
- integer(pInt), public :: &
+ !< ToDo: numerics_integrator is an array for historical reasons, only element 1 is used!
-   numerics_integrationMode   =  0_pInt                                                             !< integrationMode 1 = central solution; integrationMode 2 = perturbation, Default 0: undefined, is not read from file
+ integer(pInt), dimension(2), protected, public :: &
- integer(pInt), dimension(2) , protected, public :: &
    numerics_integrator        =  1_pInt                                                             !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states
  real(pReal), protected, public :: &
    relevantStrain             =  1.0e-7_pReal, &                                                    !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
@@ -314,9 +313,7 @@ subroutine numerics_init
        case ('atol_crystallitestress')
          aTol_crystalliteStress = IO_floatValue(line,chunkPos,2_pInt)
        case ('integrator')
-         numerics_integrator(1) = IO_intValue(line,chunkPos,2_pInt)
+         numerics_integrator = IO_intValue(line,chunkPos,2_pInt)
-       case ('integratorstiffness')
-         numerics_integrator(2) = IO_intValue(line,chunkPos,2_pInt)
        case ('usepingpong')
          usepingpong = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
        case ('timesyncing')

src/plastic_nonlocal.f90

@@ -2382,8 +2382,7 @@ use, intrinsic :: &
 use prec,     only: dNeq0, &
                     dNeq, &
                     dEq0
-use numerics, only: numerics_integrationMode, &
+use numerics, only: numerics_timeSyncing
-                    numerics_timeSyncing
 use IO,       only: IO_error
 use debug,    only: debug_level, &
                     debug_constitutive, &
@@ -2942,14 +2941,12 @@ rhoDot = rhoDotFlux &
        + rhoDotAthermalAnnihilation &
        + rhoDotThermalAnnihilation 
 
-if (numerics_integrationMode == 1_pInt) then                                                        ! save rates for output if in central integration mode
+rhoDotFluxOutput(1:ns,1:8,1_pInt,ip,el) = rhoDotFlux(1:ns,1:8)
-  rhoDotFluxOutput(1:ns,1:8,1_pInt,ip,el) = rhoDotFlux(1:ns,1:8)
+rhoDotMultiplicationOutput(1:ns,1:2,1_pInt,ip,el) = rhoDotMultiplication(1:ns,[1,3])
-  rhoDotMultiplicationOutput(1:ns,1:2,1_pInt,ip,el) = rhoDotMultiplication(1:ns,[1,3])
+rhoDotSingle2DipoleGlideOutput(1:ns,1:2,1_pInt,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10)
-  rhoDotSingle2DipoleGlideOutput(1:ns,1:2,1_pInt,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10)
+rhoDotAthermalAnnihilationOutput(1:ns,1:2,1_pInt,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10)
-  rhoDotAthermalAnnihilationOutput(1:ns,1:2,1_pInt,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10)
+rhoDotThermalAnnihilationOutput(1:ns,1:2,1_pInt,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10)
-  rhoDotThermalAnnihilationOutput(1:ns,1:2,1_pInt,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10)
+rhoDotEdgeJogsOutput(1:ns,1_pInt,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1)
-  rhoDotEdgeJogsOutput(1:ns,1_pInt,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1)
-endif
 
 
 #ifdef DEBUG