diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 4b365c68d..66ee395aa 100644 --- a/src/crystallite.f90 +++ b/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 (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 - + if (any(crystallite_todo)) call integrateState() 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, & - numerics_integrationMode + rTol_crystalliteState 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 diff --git a/src/numerics.f90 b/src/numerics.f90 index 6866a89ad..c93c4ddd7 100644 --- a/src/numerics.f90 +++ b/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 :: & - 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 :: & + !< ToDo: numerics_integrator is an array for historical reasons, only element 1 is used! + 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) - case ('integratorstiffness') - numerics_integrator(2) = IO_intValue(line,chunkPos,2_pInt) + numerics_integrator = IO_intValue(line,chunkPos,2_pInt) case ('usepingpong') usepingpong = IO_intValue(line,chunkPos,2_pInt) > 0_pInt case ('timesyncing') diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index 41666a34c..e1355da8f 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -2382,8 +2382,7 @@ use, intrinsic :: & use prec, only: dNeq0, & dNeq, & dEq0 -use numerics, only: numerics_integrationMode, & - numerics_timeSyncing +use numerics, only: 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) - 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) - 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) - rhoDotEdgeJogsOutput(1:ns,1_pInt,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1) -endif +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]) +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) +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) #ifdef DEBUG