From 73491f3be9fee6a7437306fc4b4fa2486150ff42 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Feb 2020 18:37:17 +0100 Subject: [PATCH] simplified and unified style --- src/source_damage_anisoBrittle.f90 | 97 ++++++++++------------ src/source_damage_anisoDuctile.f90 | 123 +++++++++++++--------------- src/source_damage_isoBrittle.f90 | 91 ++++++++++---------- src/source_damage_isoDuctile.f90 | 79 +++++++++--------- src/source_thermal_dissipation.f90 | 65 +++++++-------- src/source_thermal_externalheat.f90 | 66 +++++++-------- 6 files changed, 243 insertions(+), 278 deletions(-) diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index e5ed05799..2d9f41afb 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -21,14 +21,14 @@ module source_damage_anisoBrittle integer, dimension(:), allocatable :: & source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? source_damage_anisoBrittle_instance !< instance of source mechanism - + integer, dimension(:,:), allocatable :: & source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family - enum, bind(c) + enum, bind(c) enumerator :: undefined_ID, & damage_drivingforce_ID - end enum + end enum type :: tParameters !< container type for internal constitutive parameters @@ -67,68 +67,64 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoBrittle_init - integer :: Ninstance,phase,instance,source,sourceOffset - integer :: NofMyPhase,p ,i + integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p,i integer(kind(undefined_ID)) :: & outputID - character(len=pStringLen) :: & extmsg = '' character(len=pStringLen), dimension(:), allocatable :: & outputs - write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'; flush(6) + write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>'; flush(6) - Ninstance = count(phase_source == SOURCE_damage_anisoBrittle_ID) - if (Ninstance == 0) return - + Ninstance = count(phase_source == SOURCE_DAMAGE_ANISOBRITTLE_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0) - allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0) - do phase = 1, material_Nphase - source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_damage_anisoBrittle_ID) & - source_damage_anisoBrittle_offset(phase) = source - enddo - enddo - + + allocate(source_damage_anisoBrittle_offset (size(config_phase)), source=0) + allocate(source_damage_anisoBrittle_instance(size(config_phase)), source=0) + allocate(param(Ninstance)) + allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0) - allocate(param(Ninstance)) - - do p=1, size(config_phase) + do p = 1, size(config_phase) + source_damage_anisoBrittle_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ANISOBRITTLE_ID) + do source = 1, phase_Nsources(p) + if (phase_source(source,p) == SOURCE_DAMAGE_ANISOBRITTLE_ID) & + source_damage_anisoBrittle_offset(p) = source + enddo + if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISOBRITTLE_ID)) cycle + associate(prm => param(source_damage_anisoBrittle_instance(p)), & config => config_phase(p)) - + prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal) prm%N = config%getFloat('anisobrittle_ratesensitivity') prm%sdot_0 = config%getFloat('anisobrittle_sdot0') - + ! sanity checks if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol' - + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity' if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0' - + prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray) prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage)) prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage)) - + prm%cleavage_systems = lattice_SchmidMatrix_cleavage (prm%Ncleavage,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - ! expand: family => system - prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage) - prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) - - if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload' - if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement' + ! expand: family => system + prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage) + prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) + + if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload' + if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement' + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & @@ -141,7 +137,7 @@ subroutine source_damage_anisoBrittle_init do i=1, size(outputs) outputID = undefined_ID select case(outputs(i)) - + case ('anisobrittle_drivingforce') prm%outputID = [prm%outputID, damage_drivingforce_ID] @@ -150,16 +146,13 @@ subroutine source_damage_anisoBrittle_init enddo end associate - - phase = p - NofMyPhase=count(material_phaseAt==phase) * discretization_nIP - instance = source_damage_anisoBrittle_instance(phase) - sourceOffset = source_damage_anisoBrittle_offset(phase) + NofMyPhase=count(material_phaseAt==p) * discretization_nIP + instance = source_damage_anisoBrittle_instance(p) + sourceOffset = source_damage_anisoBrittle_offset(p) - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) - sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) + sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol source_damage_anisoBrittle_Ncleavage(1:size(param(instance)%Ncleavage),instance) = param(instance)%Ncleavage enddo @@ -195,9 +188,9 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) sourceOffset = source_damage_anisoBrittle_offset(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) - + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal - + index = 1 do f = 1,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family @@ -206,7 +199,7 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) - + traction_crit = param(instance)%critLoad(index)* & damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) @@ -242,12 +235,12 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph sourceOffset sourceOffset = source_damage_anisoBrittle_offset(phase) - + localphiDot = 1.0_pReal & - sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi - + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) - + end subroutine source_damage_anisoBrittle_getRateAndItsTangent @@ -257,9 +250,9 @@ end subroutine source_damage_anisoBrittle_getRateAndItsTangent subroutine source_damage_anisoBrittle_results(phase,group) integer, intent(in) :: phase - character(len=*), intent(in) :: group + character(len=*), intent(in) :: group integer :: sourceOffset, o, instance - + instance = source_damage_anisoBrittle_instance(phase) sourceOffset = source_damage_anisoBrittle_offset(phase) diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index fef897914..9e48fc343 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -13,20 +13,20 @@ module source_damage_anisoDuctile use material use config use results - + implicit none private - + integer, dimension(:), allocatable :: & source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_instance !< instance of damage source mechanism - - enum, bind(c) + + enum, bind(c) enumerator :: undefined_ID, & damage_drivingforce_ID - end enum - - + end enum + + type, private :: tParameters !< container type for internal constitutive parameters real(pReal) :: & aTol, & @@ -40,10 +40,10 @@ module source_damage_anisoDuctile integer(kind(undefined_ID)), allocatable, dimension(:) :: & outputID end type tParameters - + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) - - + + public :: & source_damage_anisoDuctile_init, & source_damage_anisoDuctile_dotState, & @@ -58,61 +58,54 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoDuctile_init - - integer :: Ninstance,phase,instance,source,sourceOffset - integer :: NofMyPhase,p ,i - + + integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p,i integer(kind(undefined_ID)) :: & outputID - character(len=pStringLen) :: & extmsg = '' character(len=pStringLen), dimension(:), allocatable :: & outputs - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'; flush(6) - - Ninstance = count(phase_source == SOURCE_damage_anisoDuctile_ID) - if (Ninstance == 0) return - + + write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init -+>>>'; flush(6) + + Ninstance = count(phase_source == SOURCE_DAMAGE_ANISODUCTILE_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_damage_anisoDuctile_offset(size(config_phase)), source=0) + + allocate(source_damage_anisoDuctile_offset (size(config_phase)), source=0) allocate(source_damage_anisoDuctile_instance(size(config_phase)), source=0) - do phase = 1, size(config_phase) - source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_damage_anisoDuctile_ID) & - source_damage_anisoDuctile_offset(phase) = source - enddo - enddo - allocate(param(Ninstance)) - - do p=1, size(config_phase) + + do p = 1, size(config_phase) + source_damage_anisoDuctile_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ANISODUCTILE_ID) + do source = 1, phase_Nsources(p) + if (phase_source(source,p) == SOURCE_DAMAGE_ANISODUCTILE_ID) & + source_damage_anisoDuctile_offset(p) = source + enddo + if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISODUCTILE_ID)) cycle + associate(prm => param(source_damage_anisoDuctile_instance(p)), & config => config_phase(p)) - + prm%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal) - + prm%N = config%getFloat('anisoductile_ratesensitivity') prm%totalNslip = sum(prm%Nslip) ! sanity checks if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol' - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity' - + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) - + prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip)) - + ! expand: family => system prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip) - + if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain' - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')') @@ -124,27 +117,25 @@ subroutine source_damage_anisoDuctile_init do i=1, size(outputs) outputID = undefined_ID select case(outputs(i)) - + case ('anisoductile_drivingforce') prm%outputID = [prm%outputID, damage_drivingforce_ID] - + end select - + enddo - + end associate - - phase = p - - NofMyPhase=count(material_phaseAt==phase) * discretization_nIP - instance = source_damage_anisoDuctile_instance(phase) - sourceOffset = source_damage_anisoDuctile_offset(phase) - - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) - sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - + + NofMyPhase=count(material_phaseAt==p) * discretization_nIP + instance = source_damage_anisoDuctile_instance(p) + sourceOffset = source_damage_anisoDuctile_offset(p) + + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) + sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol + enddo - + end subroutine source_damage_anisoDuctile_init @@ -164,22 +155,22 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) homog, damageOffset, & instance, & i - + phase = material_phaseAt(ipc,el) constituent = material_phasememberAt(ipc,ip,el) instance = source_damage_anisoDuctile_instance(phase) sourceOffset = source_damage_anisoDuctile_offset(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) - - + + do i = 1, param(instance)%totalNslip sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + & plasticState(phase)%slipRate(i,constituent)/ & - ((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i) + ((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i) enddo - + end subroutine source_damage_anisoDuctile_dotState @@ -198,14 +189,14 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph dLocalphiDot_dPhi integer :: & sourceOffset - + sourceOffset = source_damage_anisoDuctile_offset(phase) - + localphiDot = 1.0_pReal & - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi - + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) - + end subroutine source_damage_anisoDuctile_getRateAndItsTangent @@ -217,7 +208,7 @@ subroutine source_damage_anisoDuctile_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group integer :: sourceOffset, o, instance - + instance = source_damage_anisoDuctile_instance(phase) sourceOffset = source_damage_anisoDuctile_offset(phase) diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 53c0b77a7..83c4286b0 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -20,7 +20,7 @@ module source_damage_isoBrittle source_damage_isoBrittle_offset, & source_damage_isoBrittle_instance - enum, bind(c) + enum, bind(c) enumerator :: & undefined_ID, & damage_drivingforce_ID @@ -54,52 +54,47 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoBrittle_init - integer :: Ninstance,phase,instance,source,sourceOffset - integer :: NofMyPhase,p,i + integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p,i integer(kind(undefined_ID)) :: & outputID - character(len=pStringLen) :: & extmsg = '' character(len=pStringLen), dimension(:), allocatable :: & outputs - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>'; flush(6) - - Ninstance = count(phase_source == SOURCE_damage_isoBrittle_ID) - if (Ninstance == 0) return - + + write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>'; flush(6) + + Ninstance = count(phase_source == SOURCE_DAMAGE_ISOBRITTLE_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_damage_isoBrittle_offset(material_Nphase), source=0) - allocate(source_damage_isoBrittle_instance(material_Nphase), source=0) - do phase = 1, material_Nphase - source_damage_isoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoBrittle_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_damage_isoBrittle_ID) & - source_damage_isoBrittle_offset(phase) = source - enddo - enddo - + + allocate(source_damage_isoBrittle_offset (size(config_phase)), source=0) + allocate(source_damage_isoBrittle_instance(size(config_phase)), source=0) allocate(param(Ninstance)) - - do p=1, size(config_phase) + + do p = 1, size(config_phase) + source_damage_isoBrittle_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ISOBRITTLE_ID) + do source = 1, phase_Nsources(p) + if (phase_source(source,p) == SOURCE_DAMAGE_ISOBRITTLE_ID) & + source_damage_isoBrittle_offset(p) = source + enddo + if (all(phase_source(:,p) /= SOURCE_DAMAGE_ISOBRITTLE_ID)) cycle + associate(prm => param(source_damage_isoBrittle_instance(p)), & config => config_phase(p)) - + prm%aTol = config%getFloat('isobrittle_atol',defaultVal = 1.0e-3_pReal) - + prm%N = config%getFloat('isobrittle_n') prm%critStrainEnergy = config%getFloat('isobrittle_criticalstrainenergy') - + ! sanity checks if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_atol' - + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_n' if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy' - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & @@ -112,27 +107,25 @@ subroutine source_damage_isoBrittle_init do i=1, size(outputs) outputID = undefined_ID select case(outputs(i)) - + case ('isobrittle_drivingforce') prm%outputID = [prm%outputID, damage_drivingforce_ID] - + end select enddo end associate - - phase = p - - NofMyPhase = count(material_phaseAt==phase) * discretization_nIP - instance = source_damage_isoBrittle_instance(phase) - sourceOffset = source_damage_isoBrittle_offset(phase) - - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,1) - sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - + + NofMyPhase = count(material_phaseAt==p) * discretization_nIP + instance = source_damage_isoBrittle_instance(p) + sourceOffset = source_damage_isoBrittle_offset(p) + + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,1) + sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol + enddo - + end subroutine source_damage_isoBrittle_init !-------------------------------------------------------------------------------------------------- @@ -160,12 +153,12 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) instance = source_damage_isoBrittle_instance(phase) !< instance of damage_isoBrittle source sourceOffset = source_damage_isoBrittle_offset(phase) - + strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/param(instance)%critStrainEnergy ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/param(instance)%critStrainEnergy - + if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent) @@ -174,9 +167,9 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - & sourceState(phase)%p(sourceOffset)%state(1,constituent) endif - + end subroutine source_damage_isoBrittle_deltaState - + !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- @@ -195,13 +188,13 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD instance = source_damage_isoBrittle_instance(phase) sourceOffset = source_damage_isoBrittle_offset(phase) - + localphiDot = (1.0_pReal - phi)**(param(instance)%N - 1.0_pReal) - & phi*sourceState(phase)%p(sourceOffset)%state(1,constituent) dLocalphiDot_dPhi = - (param(instance)%N - 1.0_pReal)* & (1.0_pReal - phi)**max(0.0_pReal,param(instance)%N - 2.0_pReal) & - sourceState(phase)%p(sourceOffset)%state(1,constituent) - + end subroutine source_damage_isoBrittle_getRateAndItsTangent @@ -211,9 +204,9 @@ end subroutine source_damage_isoBrittle_getRateAndItsTangent subroutine source_damage_isoBrittle_results(phase,group) integer, intent(in) :: phase - character(len=*), intent(in) :: group + character(len=*), intent(in) :: group integer :: sourceOffset, o, instance - + instance = source_damage_isoBrittle_instance(phase) sourceOffset = source_damage_isoBrittle_offset(phase) diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 6ee588d0c..13d0cb3e7 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -19,7 +19,7 @@ module source_damage_isoDuctile source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_instance !< instance of damage source mechanism - enum, bind(c) + enum, bind(c) enumerator :: undefined_ID, & damage_drivingforce_ID end enum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ToDo @@ -51,52 +51,47 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoDuctile_init - integer :: Ninstance,phase,instance,source,sourceOffset - integer :: NofMyPhase,p,i + integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p,i integer(kind(undefined_ID)) :: & outputID - character(len=pStringLen) :: & extmsg = '' character(len=pStringLen), dimension(:), allocatable :: & outputs - write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init -+>>>' + write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init -+>>>'; flush(6) - Ninstance = count(phase_source == SOURCE_damage_isoDuctile_ID) - if (Ninstance == 0) return - + Ninstance = count(phase_source == SOURCE_DAMAGE_ISODUCTILE_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_damage_isoDuctile_offset(material_Nphase), source=0) - allocate(source_damage_isoDuctile_instance(material_Nphase), source=0) - do phase = 1, material_Nphase - source_damage_isoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoDuctile_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_damage_isoDuctile_ID) & - source_damage_isoDuctile_offset(phase) = source - enddo - enddo + allocate(source_damage_isoDuctile_offset (size(config_phase)), source=0) + allocate(source_damage_isoDuctile_instance(size(config_phase)), source=0) allocate(param(Ninstance)) - - do p=1, size(config_phase) + + do p = 1, size(config_phase) + source_damage_isoDuctile_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ISODUCTILE_ID) + do source = 1, phase_Nsources(p) + if (phase_source(source,p) == SOURCE_DAMAGE_ISODUCTILE_ID) & + source_damage_isoDuctile_offset(p) = source + enddo + if (all(phase_source(:,p) /= SOURCE_DAMAGE_ISODUCTILE_ID)) cycle + associate(prm => param(source_damage_isoDuctile_instance(p)), & config => config_phase(p)) - + prm%aTol = config%getFloat('isoductile_atol',defaultVal = 1.0e-3_pReal) prm%N = config%getFloat('isoductile_ratesensitivity') prm%critPlasticStrain = config%getFloat('isoductile_criticalplasticstrain') - + ! sanity checks if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isoductile_atol' - + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_ratesensitivity' if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain' - + !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') & @@ -109,7 +104,7 @@ subroutine source_damage_isoDuctile_init do i=1, size(outputs) outputID = undefined_ID select case(outputs(i)) - + case ('isoductile_drivingforce') prm%outputID = [prm%outputID, damage_drivingforce_ID] @@ -118,17 +113,17 @@ subroutine source_damage_isoDuctile_init enddo end associate - - phase = p - NofMyPhase=count(material_phaseAt==phase) * discretization_nIP - instance = source_damage_isoDuctile_instance(phase) - sourceOffset = source_damage_isoDuctile_offset(phase) - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) - sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - + + NofMyPhase=count(material_phaseAt==p) * discretization_nIP + instance = source_damage_isoDuctile_instance(p) + sourceOffset = source_damage_isoDuctile_offset(p) + + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) + sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol + enddo - + end subroutine source_damage_isoDuctile_init !-------------------------------------------------------------------------------------------------- @@ -152,11 +147,11 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & sum(plasticState(phase)%slipRate(:,constituent))/ & - ((damage(homog)%p(damageOffset))**param(instance)%N)/ & - param(instance)%critPlasticStrain + ((damage(homog)%p(damageOffset))**param(instance)%N)/ & + param(instance)%critPlasticStrain end subroutine source_damage_isoDuctile_dotState - + !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- @@ -174,12 +169,12 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD sourceOffset sourceOffset = source_damage_isoDuctile_offset(phase) - + localphiDot = 1.0_pReal & - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi - + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) - + end subroutine source_damage_isoDuctile_getRateAndItsTangent @@ -189,9 +184,9 @@ end subroutine source_damage_isoDuctile_getRateAndItsTangent subroutine source_damage_isoDuctile_results(phase,group) integer, intent(in) :: phase - character(len=*), intent(in) :: group + character(len=*), intent(in) :: group integer :: sourceOffset, o, instance - + instance = source_damage_isoDuctile_instance(phase) sourceOffset = source_damage_isoDuctile_offset(phase) diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index e13742a90..a039dbff6 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -10,26 +10,26 @@ module source_thermal_dissipation use discretization use material use config - + implicit none private integer, dimension(:), allocatable :: & source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism - + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & kappa end type tParameters - + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - - + + public :: & source_thermal_dissipation_init, & source_thermal_dissipation_getRateAndItsTangent - + contains @@ -38,45 +38,42 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_thermal_dissipation_init - - integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'; flush(6) - - - Ninstance = count(phase_source == SOURCE_thermal_dissipation_ID) - if (Ninstance == 0) return + + integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p + + write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'; flush(6) + + Ninstance = count(phase_source == SOURCE_THERMAL_DISSIPATION_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_thermal_dissipation_offset(material_Nphase), source=0) - allocate(source_thermal_dissipation_instance(material_Nphase), source=0) + + allocate(source_thermal_dissipation_offset (size(config_phase)), source=0) + allocate(source_thermal_dissipation_instance(size(config_phase)), source=0) allocate(param(Ninstance)) - - do p = 1, material_Nphase - source_thermal_dissipation_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_dissipation_ID) + + do p = 1, size(config_phase) + source_thermal_dissipation_instance(p) = count(phase_source(:,1:p) == SOURCE_THERMAL_DISSIPATION_ID) do source = 1, phase_Nsources(p) - if (phase_source(source,p) == SOURCE_thermal_dissipation_ID) & + if (phase_source(source,p) == SOURCE_THERMAL_DISSIPATION_ID) & source_thermal_dissipation_offset(p) = source - enddo - enddo - - do p=1, size(config_phase) + enddo + if (all(phase_source(:,p) /= SOURCE_THERMAL_DISSIPATION_ID)) cycle + instance = source_thermal_dissipation_instance(p) param(instance)%kappa = config_phase(p)%getFloat('dissipation_coldworkcoeff') NofMyPhase = count(material_phaseAt==p) * discretization_nIP sourceOffset = source_thermal_dissipation_offset(p) - + call material_allocateSourceState(p,sourceOffset,NofMyPhase,0,0,0) - + enddo - + end subroutine source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- -!> @brief returns dissipation rate +!> @brief Ninstances dissipation rate !-------------------------------------------------------------------------------------------------- subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase) @@ -91,12 +88,12 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar dTDOT_dT integer :: & instance - + instance = source_thermal_dissipation_instance(phase) - + TDot = param(instance)%kappa*sum(abs(Tstar*Lp)) - dTDOT_dT = 0.0_pReal - + dTDOT_dT = 0.0_pReal + end subroutine source_thermal_dissipation_getRateAndItsTangent - + end module source_thermal_dissipation diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index ccc6ebc29..182ae5f43 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -42,44 +42,40 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_init - - integer :: maxNinstance,instance,source,sourceOffset,NofMyPhase,p - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'; flush(6) - - - maxNinstance = count(phase_source == SOURCE_thermal_externalheat_ID) - if (maxNinstance == 0) return + + integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p + + write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'; flush(6) + + Ninstance = count(phase_source == SOURCE_thermal_externalheat_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(source_thermal_externalheat_offset(material_Nphase), source=0) - allocate(source_thermal_externalheat_instance(material_Nphase), source=0) - - do p = 1, material_Nphase + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(source_thermal_externalheat_offset (size(config_phase)), source=0) + allocate(source_thermal_externalheat_instance(size(config_phase)), source=0) + allocate(param(Ninstance)) + + do p = 1, size(config_phase) source_thermal_externalheat_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_externalheat_ID) do source = 1, phase_Nsources(p) if (phase_source(source,p) == SOURCE_thermal_externalheat_ID) & source_thermal_externalheat_offset(p) = source - enddo - enddo - - allocate(param(maxNinstance)) - - do p=1, size(config_phase) + enddo + if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle + instance = source_thermal_externalheat_instance(p) sourceOffset = source_thermal_externalheat_offset(p) NofMyPhase = count(material_phaseAt==p) * discretization_nIP - + param(instance)%time = config_phase(p)%getFloats('externalheat_time') param(instance)%nIntervals = size(param(instance)%time) - 1 - - + + param(instance)%heat_rate = config_phase(p)%getFloats('externalheat_rate',requiredSize = size(param(instance)%time)) - + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) - + enddo end subroutine source_thermal_externalheat_init @@ -90,25 +86,25 @@ end subroutine source_thermal_externalheat_init !> @details state only contains current time to linearly interpolate given heat powers !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_dotState(phase, of) - + integer, intent(in) :: & phase, & of integer :: & sourceOffset - + sourceOffset = source_thermal_externalheat_offset(phase) - + sourceState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- -!> @brief returns local heat generation rate +!> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) - + integer, intent(in) :: & phase, & of @@ -118,11 +114,11 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas integer :: & instance, sourceOffset, interval real(pReal) :: & - frac_time - + frac_time + instance = source_thermal_externalheat_instance(phase) sourceOffset = source_thermal_externalheat_offset(phase) - + do interval = 1, param(instance)%nIntervals ! scan through all rate segments frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - param(instance)%time(interval)) & / (param(instance)%time(interval+1) - param(instance)%time(interval)) ! fractional time within segment @@ -134,7 +130,7 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas ! ...or extrapolate if outside of bounds enddo dTDot_dT = 0.0 - + end subroutine source_thermal_externalheat_getRateAndItsTangent - + end module source_thermal_externalheat