diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 14b624b0a..eb10e8550 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -23,6 +23,7 @@ module crystallite use discretization use lattice use results + use YAML_types implicit none private @@ -81,8 +82,9 @@ module crystallite integer :: & iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp nState, & !< state loop limit - nStress, & !< stress loop limit - integrator !< integration scheme (ToDo: better use a string) + nStress !< stress loop limit + character(len=pStringLen) :: & + integrator !< integrator scheme real(pReal) :: & subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback subStepSizeCryst, & !< size of first substep when cutback @@ -126,7 +128,9 @@ subroutine crystallite_init iMax, & !< maximum number of integration points eMax, & !< maximum number of elements myNcomponents !< number of components at current IP - + + class(tNode) , pointer :: & + num_crystallite write(6,'(/,a)') ' <<<+- crystallite init -+>>>' cMax = homogenization_maxNgrains @@ -160,23 +164,20 @@ subroutine crystallite_init allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) - num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal) - num%subStepSizeCryst = config_numerics%getFloat('substepsizecryst', defaultVal=0.25_pReal) - num%stepIncreaseCryst = config_numerics%getFloat('stepincreasecryst', defaultVal=1.5_pReal) - - num%subStepSizeLp = config_numerics%getFloat('substepsizelp', defaultVal=0.5_pReal) - num%subStepSizeLi = config_numerics%getFloat('substepsizeli', defaultVal=0.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) - - num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1) - - num%integrator = config_numerics%getInt ('integrator', defaultVal=1) - - num%nState = config_numerics%getInt ('nstate', defaultVal=20) - num%nStress = config_numerics%getInt ('nstress', defaultVal=40) + num_crystallite => numerics_root%get('crystallite',defaultVal=emptyDict) + + num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal) + num%subStepSizeCryst = num_crystallite%get_asFloat ('subStepSize', defaultVal=0.25_pReal) + num%stepIncreaseCryst = num_crystallite%get_asFloat ('stepIncrease', defaultVal=1.5_pReal) + num%subStepSizeLp = num_crystallite%get_asFloat ('subStepSizeLp', defaultVal=0.5_pReal) + num%subStepSizeLi = num_crystallite%get_asFloat ('subStepSizeLi', defaultVal=0.5_pReal) + num%rtol_crystalliteState = num_crystallite%get_asFloat ('rtol_State', defaultVal=1.0e-6_pReal) + num%rtol_crystalliteStress = num_crystallite%get_asFloat ('rtol_Stress', defaultVal=1.0e-6_pReal) + num%atol_crystalliteStress = num_crystallite%get_asFloat ('atol_Stress', defaultVal=1.0e-8_pReal) + num%iJacoLpresiduum = num_crystallite%get_asInt ('iJacoLpresiduum', defaultVal=1) + num%integrator = num_crystallite%get_asString('integrator', defaultVal='FPI') + num%nState = num_crystallite%get_asInt ('nState', defaultVal=20) + num%nStress = num_crystallite%get_asInt ('nStress', defaultVal=40) 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') @@ -191,24 +192,23 @@ subroutine crystallite_init if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') - if(num%integrator < 1 .or. num%integrator > 5) & - call IO_error(301,ext_msg='integrator') - if(num%nState < 1) call IO_error(301,ext_msg='nState') if(num%nStress< 1) call IO_error(301,ext_msg='nStress') - select case(num%integrator) - case(1) + select case(trim(num%integrator)) + case('FPI') integrateState => integrateStateFPI - case(2) + case('Euler') integrateState => integrateStateEuler - case(3) + case('AdaptiveEuler') integrateState => integrateStateAdaptiveEuler - case(4) + case('RK4') integrateState => integrateStateRK4 - case(5) + case('RKCK45') integrateState => integrateStateRKCK45 + case default + call IO_error(301,ext_msg='integrator') end select allocate(output_constituent(size(config_phase))) @@ -245,7 +245,6 @@ subroutine crystallite_init enddo !$OMP END PARALLEL DO - !if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFi0 = crystallite_Fi0 @@ -276,6 +275,7 @@ subroutine crystallite_init write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax flush(6) endif + #endif end subroutine crystallite_init diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 1c2f51056..93117ef26 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -92,9 +92,9 @@ function damage_local_updateState(subdt, ip, el) phi = max(residualStiffness,min(1.0_pReal,phi + subdt*phiDot)) damage_local_updateState = [ abs(phi - damageState(homog)%state(1,offset)) & - <= err_damage_tolAbs & + <= 1.0e-2_pReal & .or. abs(phi - damageState(homog)%state(1,offset)) & - <= err_damage_tolRel*abs(damageState(homog)%state(1,offset)), & + <= 1.0e-6_pReal*abs(damageState(homog)%state(1,offset)), & .true.] damageState(homog)%state(1,offset) = phi diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index eda60d070..41ed79f40 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -165,7 +165,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) phi_stagInc = phi_current - solution%stagConverged = stagNorm < max(err_damage_tolAbs, err_damage_tolRel*solnNorm) + solution%stagConverged = stagNorm < max(1.0e-2_pReal, 1.0e-6_pReal*solnNorm) !-------------------------------------------------------------------------------------------------- ! updating damage state diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 1c8314369..40e27ce34 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -162,7 +162,7 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) T_stagInc = T_current - solution%stagConverged = stagNorm < max(err_thermal_tolAbs, err_thermal_tolRel*solnNorm) + solution%stagConverged = stagNorm < max(1.0e-2_pReal, 1.0e-6_pReal*solnNorm) !-------------------------------------------------------------------------------------------------- ! updating thermal state diff --git a/src/numerics.f90 b/src/numerics.f90 index 77c3397de..c677ef888 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -6,6 +6,8 @@ module numerics use prec use IO + use YAML_types + use YAML_parse #ifdef PETSc #include @@ -16,6 +18,8 @@ module numerics implicit none private + class(tNode), pointer, public :: & + numerics_root integer, protected, public :: & iJacoStiffness = 1, & !< frequency of stiffness update randomSeed = 0, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed @@ -33,11 +37,7 @@ module numerics ! field parameters: real(pReal), protected, public :: & err_struct_tolAbs = 1.0e-10_pReal, & !< absolute tolerance for mechanical equilibrium - err_struct_tolRel = 1.0e-4_pReal, & !< relative tolerance for mechanical equilibrium - err_thermal_tolAbs = 1.0e-2_pReal, & !< absolute tolerance for thermal equilibrium - err_thermal_tolRel = 1.0e-6_pReal, & !< relative tolerance for thermal equilibrium - err_damage_tolAbs = 1.0e-2_pReal, & !< absolute tolerance for damage evolution - err_damage_tolRel = 1.0e-6_pReal !< relative tolerance for damage evolution + err_struct_tolRel = 1.0e-4_pReal !< relative tolerance for mechanical equilibrium integer, protected, public :: & itmax = 250, & !< maximum number of iterations itmin = 1, & !< minimum number of iterations @@ -84,11 +84,14 @@ contains subroutine numerics_init !$ integer :: gotDAMASK_NUM_THREADS = 1 integer :: i,j, ierr - integer, allocatable, dimension(:) :: chunkPos - character(len=pStringLen), dimension(:), allocatable :: fileContent - character(len=pStringLen) :: & - tag ,& - line + character(len=:), allocatable :: & + numerics_input, & + numerics_inFlow, & + key + class (tNode), pointer :: & + num_grid, & + num_mesh, & + num_generic logical :: fexist !$ character(len=6) DAMASK_NumThreadsString ! environment variable DAMASK_NUM_THREADS @@ -108,111 +111,104 @@ subroutine numerics_init !$ endif !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution - inquire(file='numerics.config', exist=fexist) - + inquire(file='numerics.yaml', exist=fexist) + fileExists: if (fexist) then write(6,'(a,/)') ' using values from config file' flush(6) - fileContent = IO_read_ASCII('numerics.config') - do j=1, size(fileContent) - + numerics_input = IO_read('numerics.yaml') + numerics_inFlow = to_flow(numerics_input) + numerics_root => parse_flow(numerics_inFlow,defaultVal=emptyDict) !-------------------------------------------------------------------------------------------------- -! read variables from config file and overwrite default parameters if keyword is present - line = fileContent(j) - do i=1,len(line) - if(line(i:i) == '=') line(i:i) = ' ' ! also allow keyword = value version - enddo - if (IO_isBlank(line)) cycle ! skip empty lines - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1)) ! extract key - - select case(tag) +! spectral parameters + num_grid => numerics_root%get('grid',defaultVal=emptyDict) + do i=1,num_grid%length + key = num_grid%getKey(i) + select case(key) +#ifdef Grid + case ('err_div_tolabs') + err_div_tolAbs = num_grid%get_asFloat(key) + case ('err_div_tolrel') + err_div_tolRel = num_grid%get_asFloat(key) + case ('err_stress_tolrel') + err_stress_tolrel = num_grid%get_asFloat(key) + case ('err_stress_tolabs') + err_stress_tolabs = num_grid%get_asFloat(key) + case ('err_curl_tolabs') + err_curl_tolAbs = num_grid%get_asFloat(key) + case ('err_curl_tolrel') + err_curl_tolRel = num_grid%get_asFloat(key) + case ('polaralpha') + polarAlpha = num_grid%get_asFloat(key) + case ('polarbeta') + polarBeta = num_grid%get_asFloat(key) +#endif + case ('itmax') + itmax = num_grid%get_asInt(key) + case ('itmin') + itmin = num_grid%get_asInt(key) + case ('maxCutBack') + maxCutBack = num_grid%get_asInt(key) + case ('maxStaggeredIter') + stagItMax = num_grid%get_asInt(key) +#ifdef PETSC + case ('petsc_options') + petsc_options = num_grid%get_asString(key) +#endif + endselect + enddo + + num_generic => numerics_root%get('generic',defaultVal=emptyDict) + do i=1,num_generic%length + key = num_generic%getKey(i) + select case(key) case ('defgradtolerance') - defgradTolerance = IO_floatValue(line,chunkPos,2) + defgradTolerance = num_generic%get_asFloat(key) case ('ijacostiffness') - iJacoStiffness = IO_intValue(line,chunkPos,2) + iJacoStiffness = num_generic%get_asInt(key) case ('unitlength') - numerics_unitlength = IO_floatValue(line,chunkPos,2) + numerics_unitlength = num_generic%get_asFloat(key) !-------------------------------------------------------------------------------------------------- ! random seeding parameter - case ('random_seed','fixed_seed') - randomSeed = IO_intValue(line,chunkPos,2) + case ('fixed_seed', 'random_seed') + randomSeed = num_generic%get_asInt(key) !-------------------------------------------------------------------------------------------------- ! gradient parameter - case ('charlength') - charLength = IO_floatValue(line,chunkPos,2) - case ('residualstiffness') - residualStiffness = IO_floatValue(line,chunkPos,2) - + case ('charLength') + charLength = num_generic%get_asFloat(key) + case ('residualStiffness') + residualStiffness = num_generic%get_asFloat(key) !-------------------------------------------------------------------------------------------------- ! field parameters case ('err_struct_tolabs') - err_struct_tolAbs = IO_floatValue(line,chunkPos,2) + err_struct_tolAbs = num_generic%get_asFloat(key) case ('err_struct_tolrel') - err_struct_tolRel = IO_floatValue(line,chunkPos,2) - case ('err_thermal_tolabs') - err_thermal_tolabs = IO_floatValue(line,chunkPos,2) - case ('err_thermal_tolrel') - err_thermal_tolrel = IO_floatValue(line,chunkPos,2) - case ('err_damage_tolabs') - err_damage_tolabs = IO_floatValue(line,chunkPos,2) - case ('err_damage_tolrel') - err_damage_tolrel = IO_floatValue(line,chunkPos,2) - case ('itmax') - itmax = IO_intValue(line,chunkPos,2) - case ('itmin') - itmin = IO_intValue(line,chunkPos,2) - case ('maxcutback') - maxCutBack = IO_intValue(line,chunkPos,2) - case ('maxstaggerediter') - stagItMax = IO_intValue(line,chunkPos,2) + err_struct_tolRel = num_generic%get_asFloat(key) + endselect + enddo -#ifdef PETSC - case ('petsc_options') - petsc_options = trim(line(chunkPos(4):)) -#endif - -!-------------------------------------------------------------------------------------------------- -! spectral parameters -#ifdef Grid - case ('err_div_tolabs') - err_div_tolAbs = IO_floatValue(line,chunkPos,2) - case ('err_div_tolrel') - err_div_tolRel = IO_floatValue(line,chunkPos,2) - case ('err_stress_tolrel') - err_stress_tolrel = IO_floatValue(line,chunkPos,2) - case ('err_stress_tolabs') - err_stress_tolabs = IO_floatValue(line,chunkPos,2) - case ('err_curl_tolabs') - err_curl_tolAbs = IO_floatValue(line,chunkPos,2) - case ('err_curl_tolrel') - err_curl_tolRel = IO_floatValue(line,chunkPos,2) - case ('polaralpha') - polarAlpha = IO_floatValue(line,chunkPos,2) - case ('polarbeta') - polarBeta = IO_floatValue(line,chunkPos,2) -#endif - -!-------------------------------------------------------------------------------------------------- -! Mesh parameters #ifdef Mesh + num_grid => numerics_root%get('mesh',defaultVal=emptyDict) + do i=1,num_grid%length + key = num_grid%getKey(i) + select case(key) case ('integrationorder') - integrationorder = IO_intValue(line,chunkPos,2) + integrationorder = num_generic%get_asInt(key) case ('structorder') - structorder = IO_intValue(line,chunkPos,2) + structorder = num_generic%get_asInt(key) case ('bbarstabilisation') - BBarStabilisation = IO_intValue(line,chunkPos,2) > 0 -#endif + BBarStabilisation = num_generic%get_asInt(key) > 0 end select enddo +#endif + else fileExists write(6,'(a,/)') ' using standard values' flush(6) endif fileExists - !-------------------------------------------------------------------------------------------------- ! writing parameters to output write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance @@ -242,10 +238,6 @@ subroutine numerics_init write(6,'(a24,1x,i8)') ' maxStaggeredIter: ',stagItMax write(6,'(a24,1x,es8.1)') ' err_struct_tolAbs: ',err_struct_tolAbs write(6,'(a24,1x,es8.1)') ' err_struct_tolRel: ',err_struct_tolRel - write(6,'(a24,1x,es8.1)') ' err_thermal_tolabs: ',err_thermal_tolabs - write(6,'(a24,1x,es8.1)') ' err_thermal_tolrel: ',err_thermal_tolrel - write(6,'(a24,1x,es8.1)') ' err_damage_tolabs: ',err_damage_tolabs - write(6,'(a24,1x,es8.1)') ' err_damage_tolrel: ',err_damage_tolrel !-------------------------------------------------------------------------------------------------- ! spectral parameters @@ -284,10 +276,6 @@ subroutine numerics_init if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (err_struct_tolRel <= 0.0_pReal) call IO_error(301,ext_msg='err_struct_tolRel') if (err_struct_tolAbs <= 0.0_pReal) call IO_error(301,ext_msg='err_struct_tolAbs') - if (err_thermal_tolabs <= 0.0_pReal) call IO_error(301,ext_msg='err_thermal_tolabs') - if (err_thermal_tolrel <= 0.0_pReal) call IO_error(301,ext_msg='err_thermal_tolrel') - if (err_damage_tolabs <= 0.0_pReal) call IO_error(301,ext_msg='err_damage_tolabs') - if (err_damage_tolrel <= 0.0_pReal) call IO_error(301,ext_msg='err_damage_tolrel') #ifdef Grid if (err_stress_tolrel <= 0.0_pReal) call IO_error(301,ext_msg='err_stress_tolRel') if (err_stress_tolabs <= 0.0_pReal) call IO_error(301,ext_msg='err_stress_tolAbs') diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 2bef1c3ee..a11e7b0a1 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -101,9 +101,9 @@ function thermal_adiabatic_updateState(subdt, ip, el) T = T + subdt*Tdot/(thermal_adiabatic_getSpecificHeat(ip,el)*thermal_adiabatic_getMassDensity(ip,el)) thermal_adiabatic_updateState = [ abs(T - thermalState(homog)%state(1,offset)) & - <= err_thermal_tolAbs & + <= 1.0e-2_pReal & .or. abs(T - thermalState(homog)%state(1,offset)) & - <= err_thermal_tolRel*abs(thermalState(homog)%state(1,offset)), & + <= 1.0e-6_pReal*abs(thermalState(homog)%state(1,offset)), & .true.] temperature (homog)%p(thermalMapping(homog)%p(ip,el)) = T