From 67eb39255ac29a30b02bcfdbb58ef59fbcf04c8f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Apr 2019 07:24:04 +0200 Subject: [PATCH] keep connected data together - avoids dependencies - easier to read and modify --- src/crystallite.f90 | 73 ++++++++++++++++++++++++++++++--------------- src/numerics.f90 | 32 -------------------- 2 files changed, 49 insertions(+), 56 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 69c7839c7..afaffb07a 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -111,6 +111,20 @@ module crystallite end type tOutput type(tOutput), allocatable, dimension(:), private :: output_constituent + type, private :: tNumerics + real(pReal) :: & + subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback + subStepSizeCryst, & !< size of first substep when cutback + subStepSizeLp, & !< size of first substep when cutback in Lp calculation + subStepSizeLi, & !< size of first substep when cutback in Li calculation + stepIncreaseCryst, & !< increase of next substep size when previous substep converged + rTol_crystalliteState, & !< relative tolerance in state loop + rTol_crystalliteStress, & !< relative tolerance in stress loop + aTol_crystalliteStress !< absolute tolerance in stress loop + end type tNumerics + + type(tNumerics) :: num ! numerics parameters. Better name? + procedure(), pointer :: integrateState public :: & @@ -165,6 +179,7 @@ subroutine crystallite_init use config, only: & config_deallocate, & config_crystallite, & + config_numerics, & config_phase, & crystallite_name use constitutive, only: & @@ -242,6 +257,24 @@ subroutine crystallite_init allocate(crystallite_sizePostResults(size(config_crystallite)),source=0) allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & size(config_crystallite)), source=0) + + num%subStepMinCryst = config_numerics%getFloat('subStepMinCryst', defaultVal=1.0e-3_pReal) + num%subStepSizeCryst = config_numerics%getFloat('subStepSizeCryst', defaultVal=0.25_pReal) + num%subStepSizeLp = config_numerics%getFloat('subStepSizeLp', defaultVal=0.5_pReal) + num%subStepSizeLi = config_numerics%getFloat('subStepSizeLi', defaultVal=0.5_pReal) + num%stepIncreaseCryst = config_numerics%getFloat('stepIncreaseCryst', defaultVal=1.5_pReal) + num%rTol_crystalliteState = config_numerics%getFloat('rTol_crystalliteState', defaultVal=1.0e-6_pReal) + num%rTol_crystalliteStress = config_numerics%getFloat('rTol_crystalliteStress',defaultVal=1.0e-6_pReal) + num%aTol_crystalliteStress = config_numerics%getFloat('aTol_crystalliteStress',defaultVal=1.0e-8_pReal) + + if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst') + if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst') + if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst') + if(num%subStepSizeLp <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLp') + if(num%subStepSizeLi <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLi') + if(num%rTol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteState') + if(num%rTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteStress') + if(num%aTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='aTol_crystalliteStress') select case(numerics_integrator) case(1) @@ -433,10 +466,6 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) use prec, only: & tol_math_check, & dNeq0 - use numerics, only: & - subStepMinCryst, & - subStepSizeCryst, & - stepIncreaseCryst #ifdef DEBUG use debug, only: & debug_level, & @@ -519,7 +548,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) crystallite_subS0(1:3,1:3,c,i,e) = crystallite_partionedS0(1:3,1:3,c,i,e) crystallite_subFrac(c,i,e) = 0.0_pReal - crystallite_subStep(c,i,e) = 1.0_pReal/subStepSizeCryst + crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst crystallite_todo(c,i,e) = .true. crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst endif homogenizationRequestsCalculation @@ -554,7 +583,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) formerSubStep = crystallite_subStep(c,i,e) crystallite_subFrac(c,i,e) = crystallite_subFrac(c,i,e) + crystallite_subStep(c,i,e) crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & - stepIncreaseCryst * crystallite_subStep(c,i,e)) + num%stepIncreaseCryst * crystallite_subStep(c,i,e)) crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? if (crystallite_todo(c,i,e)) then @@ -584,7 +613,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) else - crystallite_subStep(c,i,e) = subStepSizeCryst * crystallite_subStep(c,i,e) + crystallite_subStep(c,i,e) = num%subStepSizeCryst * crystallite_subStep(c,i,e) crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp (1:3,1:3,c,i,e)) crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) @@ -602,7 +631,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) enddo ! cant restore dotState here, since not yet calculated in first cutback after initialization - crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst ! still on track or already done (beyond repair) + crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & @@ -652,7 +681,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) !-------------------------------------------------------------------------------------------------- ! integrate --- requires fully defined state array (basic + dependent state) if (any(crystallite_todo)) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation - where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & ! do not try non-converged but fully cutbacked any further + where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation @@ -1238,11 +1267,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) use prec, only: tol_math_check, & dEq0 use numerics, only: nStress, & - aTol_crystalliteStress, & - rTol_crystalliteStress, & - iJacoLpresiduum, & - subStepSizeLp, & - subStepSizeLi + iJacoLpresiduum #ifdef DEBUG use debug, only: debug_level, & debug_e, & @@ -1442,8 +1467,8 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif !* update current residuum and check for convergence of loop - aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error - aTol_crystalliteStress) ! minimum lower cutoff + aTolLp = max(num%rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error + num%aTol_crystalliteStress) ! minimum lower cutoff residuumLp = Lpguess - Lp_constitutive if (any(IEEE_is_NaN(residuumLp))) then @@ -1463,7 +1488,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) Lpguess_old = Lpguess steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) else ! not converged and residuum not improved... - steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction + steplengthLp = num%subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction Lpguess = Lpguess_old + steplengthLp * deltaLp #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & @@ -1541,8 +1566,8 @@ logical function integrateStress(ipc,ip,el,timeFraction) #endif !* update current residuum and check for convergence of loop - aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error - aTol_crystalliteStress) ! minimum lower cutoff + aTolLi = max(num%rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error + num%aTol_crystalliteStress) ! minimum lower cutoff residuumLi = Liguess - Li_constitutive if (any(IEEE_is_NaN(residuumLi))) then ! NaN in residuum... #ifdef DEBUG @@ -1561,7 +1586,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) Liguess_old = Liguess steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) else ! not converged and residuum not improved... - steplengthLi = subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction + steplengthLi = num%subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction Liguess = Liguess_old + steplengthLi * deltaLi cycle LiLoop endif @@ -2295,14 +2320,14 @@ end subroutine setConvergenceFlag !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- logical pure function converged(residuum,state,aTol) - use prec, only: & - dEq0 - use numerics, only: & - rTol => rTol_crystalliteState implicit none real(pReal), intent(in), dimension(:) ::& residuum, state, aTol + real(pReal) :: & + rTol + + rTol = num%rTol_crystalliteState converged = all(abs(residuum) <= max(aTol, rTol*abs(state))) diff --git a/src/numerics.f90 b/src/numerics.f90 index 2a9920a6d..a5368a5de 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -26,17 +26,9 @@ module numerics DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive real(pReal), protected, public :: & defgradTolerance = 1.0e-7_pReal, & !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1) - subStepMinCryst = 1.0e-3_pReal, & !< minimum (relative) size of sub-step allowed during cutback in crystallite subStepMinHomog = 1.0e-3_pReal, & !< minimum (relative) size of sub-step allowed during cutback in homogenization - subStepSizeCryst = 0.25_pReal, & !< size of first substep when cutback in crystallite subStepSizeHomog = 0.25_pReal, & !< size of first substep when cutback in homogenization - subStepSizeLp = 0.5_pReal, & !< size of first substep when cutback in Lp calculation - subStepSizeLi = 0.5_pReal, & !< size of first substep when cutback in Li calculation - stepIncreaseCryst = 1.5_pReal, & !< increase of next substep size when previous substep converged in crystallite stepIncreaseHomog = 1.5_pReal, & !< increase of next substep size when previous substep converged in homogenization - rTol_crystalliteState = 1.0e-6_pReal, & !< relative tolerance in crystallite state loop - rTol_crystalliteStress = 1.0e-6_pReal, & !< relative tolerance in crystallite stress loop - aTol_crystalliteStress = 1.0e-8_pReal, & !< absolute tolerance in crystallite stress loop, Default 1.0e-8: residuum is in Lp and hence strain is on this order numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit absTol_RGC = 1.0e+4_pReal, & !< absolute tolerance of RGC residuum relTol_RGC = 1.0e-3_pReal, & !< relative tolerance of RGC residuum @@ -205,15 +197,10 @@ subroutine numerics_init case ('nstress') nStress = IO_intValue(line,chunkPos,2_pInt) case ('substepmincryst') - subStepMinCryst = IO_floatValue(line,chunkPos,2_pInt) case ('substepsizecryst') - subStepSizeCryst = IO_floatValue(line,chunkPos,2_pInt) case ('stepincreasecryst') - stepIncreaseCryst = IO_floatValue(line,chunkPos,2_pInt) case ('substepsizelp') - subStepSizeLp = IO_floatValue(line,chunkPos,2_pInt) case ('substepsizeli') - subStepSizeLi = IO_floatValue(line,chunkPos,2_pInt) case ('substepminhomog') subStepMinHomog = IO_floatValue(line,chunkPos,2_pInt) case ('substepsizehomog') @@ -221,11 +208,8 @@ subroutine numerics_init case ('stepincreasehomog') stepIncreaseHomog = IO_floatValue(line,chunkPos,2_pInt) case ('rtol_crystallitestate') - rTol_crystalliteState = IO_floatValue(line,chunkPos,2_pInt) case ('rtol_crystallitestress') - rTol_crystalliteStress = IO_floatValue(line,chunkPos,2_pInt) case ('atol_crystallitestress') - aTol_crystalliteStress = IO_floatValue(line,chunkPos,2_pInt) case ('integrator') numerics_integrator = IO_intValue(line,chunkPos,2_pInt) case ('usepingpong') @@ -350,16 +334,8 @@ subroutine numerics_init write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness write(6,'(a24,1x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum write(6,'(a24,1x,i8)') ' nCryst: ',nCryst - write(6,'(a24,1x,es8.1)') ' subStepMinCryst: ',subStepMinCryst - write(6,'(a24,1x,es8.1)') ' subStepSizeCryst: ',subStepSizeCryst - write(6,'(a24,1x,es8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst - write(6,'(a24,1x,es8.1)') ' subStepSizeLp: ',subStepSizeLp - write(6,'(a24,1x,es8.1)') ' subStepSizeLi: ',subStepSizeLi write(6,'(a24,1x,i8)') ' nState: ',nState write(6,'(a24,1x,i8)') ' nStress: ',nStress - write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState - write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress - write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,1x,i8)') ' integrator: ',numerics_integrator write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength @@ -446,17 +422,9 @@ subroutine numerics_init if (nCryst < 1_pInt) call IO_error(301_pInt,ext_msg='nCryst') if (nState < 1_pInt) call IO_error(301_pInt,ext_msg='nState') if (nStress < 1_pInt) call IO_error(301_pInt,ext_msg='nStress') - if (subStepMinCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinCryst') - if (subStepSizeCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeCryst') - if (stepIncreaseCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseCryst') - if (subStepSizeLp <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLp') - if (subStepSizeLi <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeLi') if (subStepMinHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinHomog') if (subStepSizeHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeHomog') if (stepIncreaseHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseHomog') - if (rTol_crystalliteState <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteState') - if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteStress') - if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='aTol_crystalliteStress') if (numerics_integrator <= 0_pInt .or. numerics_integrator >= 6_pInt) & call IO_error(301_pInt,ext_msg='integrator') if (numerics_unitlength <= 0.0_pReal) call IO_error(301_pInt,ext_msg='unitlength')