From 991d0fe020c3393b30062ceaa8dd5fd11510b61d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 28 Feb 2020 10:25:07 +0100 Subject: [PATCH] polishing/unifying --- src/kinematics_cleavage_opening.f90 | 35 +++++------ src/kinematics_slipplane_opening.f90 | 86 +++++++++++++++------------- src/kinematics_thermal_expansion.f90 | 62 ++++++++++---------- src/source_damage_anisoBrittle.f90 | 23 +++----- src/source_damage_anisoDuctile.f90 | 22 +++---- src/source_damage_isoBrittle.f90 | 18 +++--- src/source_damage_isoDuctile.f90 | 7 +-- src/source_thermal_dissipation.f90 | 10 ++-- src/source_thermal_externalheat.f90 | 6 +- 9 files changed, 131 insertions(+), 138 deletions(-) diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index c088a6c00..e63ea0d7d 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -63,40 +63,41 @@ subroutine kinematics_cleavage_opening_init integer, allocatable, dimension(:) :: tempInt real(pReal), allocatable, dimension(:) :: tempFloat - integer :: maxNinstance,p,instance + integer :: Ninstance,p,instance write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>'; flush(6) - maxNinstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID) - if (maxNinstance == 0) return - + Ninstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0) do p = 1, size(config_phase) kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct? enddo - allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) - allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) - allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0) - allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0) - allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal) - allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal) + allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,Ninstance), source=0.0_pReal) + allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,Ninstance), source=0.0_pReal) + allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0) + allocate(kinematics_cleavage_opening_totalNcleavage(Ninstance), source=0) + allocate(kinematics_cleavage_opening_sdot_0(Ninstance), source=0.0_pReal) + allocate(kinematics_cleavage_opening_N(Ninstance), source=0.0_pReal) do p = 1, size(config_phase) if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle + associate(config => config_phase(p)) + instance = kinematics_cleavage_opening_instance(p) - kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0') - kinematics_cleavage_opening_N(instance) = config_phase(p)%getFloat('anisobrittle_ratesensitivity') - tempInt = config_phase(p)%getInts('ncleavage') + + kinematics_cleavage_opening_sdot_0(instance) = config%getFloat('anisobrittle_sdot0') + kinematics_cleavage_opening_N(instance) = config%getFloat('anisobrittle_ratesensitivity') + tempInt = config%getInts('ncleavage') kinematics_cleavage_opening_Ncleavage(1:size(tempInt),instance) = tempInt - tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt)) + tempFloat = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt)) kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat - tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt)) + tempFloat = config%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt)) kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = & @@ -111,6 +112,8 @@ subroutine kinematics_cleavage_opening_init call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')') if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) & call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')') + + end associate enddo end subroutine kinematics_cleavage_opening_init diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 2c94448bd..48aef3871 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -12,12 +12,12 @@ module kinematics_slipplane_opening use math use lattice use material - + implicit none private - + integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance - + type :: tParameters !< container type for internal constitutive parameters integer :: & totalNslip @@ -33,9 +33,9 @@ module kinematics_slipplane_opening slip_normal, & slip_transverse end type tParameters - + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - + public :: & kinematics_slipplane_opening_init, & kinematics_slipplane_opening_LiAndItsTangent @@ -49,43 +49,39 @@ contains !-------------------------------------------------------------------------------------------------- subroutine kinematics_slipplane_opening_init - integer :: maxNinstance,p,instance + integer :: Ninstance,p,instance write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>'; flush(6) - maxNinstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID) - if (maxNinstance == 0) return - + Ninstance = count(phase_kinematics == KINEMATICS_slipplane_opening_ID) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0) do p = 1, size(config_phase) kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct? enddo - - allocate(param(maxNinstance)) - + + allocate(param(Ninstance)) + do p = 1, size(config_phase) if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle associate(prm => param(kinematics_slipplane_opening_instance(p)), & config => config_phase(p)) - instance = kinematics_slipplane_opening_instance(p) - prm%sdot0 = config_phase(p)%getFloat('anisoductile_sdot0') - prm%n = config_phase(p)%getFloat('anisoductile_ratesensitivity') - - prm%Nslip = config%getInts('nslip') - prm%critLoad = config_phase(p)%getFloats('anisoductile_criticalload',requiredSize=size(prm%Nslip )) + prm%sdot0 = config%getFloat('anisoductile_sdot0') + prm%n = config%getFloat('anisoductile_ratesensitivity') + prm%Nslip = config%getInts('nslip') + prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + + prm%critLoad = config%getFloats('anisoductile_criticalload',requiredSize=size(prm%Nslip)) prm%critLoad = math_expand(prm%critLoad, prm%Nslip) - - prm%slip_direction = lattice_slip_direction (prm%Nslip,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%slip_normal = lattice_slip_normal (prm%Nslip,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%slip_transverse = lattice_slip_transverse(prm%Nslip,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) ! if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) & ! call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')') @@ -93,14 +89,14 @@ subroutine kinematics_slipplane_opening_init ! call IO_error(211,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')') ! if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) & ! call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')') - + end associate enddo - + end subroutine kinematics_slipplane_opening_init !-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient +!> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) @@ -123,12 +119,12 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - + phase = material_phaseAt(ipc,el) instance = kinematics_slipplane_opening_instance(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) - + associate(prm => param(instance)) Ld = 0.0_pReal dLd_dTstar = 0.0_pReal @@ -137,12 +133,12 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, projection_d = math_outer(prm%slip_direction(1:3,i), prm%slip_normal(1:3,i)) projection_t = math_outer(prm%slip_transverse(1:3,i),prm%slip_normal(1:3,i)) projection_n = math_outer(prm%slip_normal(1:3,i), prm%slip_normal(1:3,i)) - + traction_d = math_mul33xx33(S,projection_d) traction_t = math_mul33xx33(S,projection_t) traction_n = math_mul33xx33(S,projection_n) - - traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage + + traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage udotd = sign(1.0_pReal,traction_d)* prm%sdot0* ( abs(traction_d)/traction_crit & - abs(traction_d)/prm%critLoad(i))**prm%n @@ -151,9 +147,21 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, udotn = prm%sdot0* ( max(0.0_pReal,traction_n)/traction_crit & - max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n - dudotd_dt = udotd*prm%n/traction_d - dudott_dt = udott*prm%n/traction_t - dudotn_dt = udotn*prm%n/traction_n + if (dNeq0(traction_d)) then + dudotd_dt = udotd*prm%n/traction_d + else + dudotd_dt = 0.0_pReal + endif + if (dNeq0(traction_t)) then + dudott_dt = udott*prm%n/traction_t + else + dudott_dt = 0.0_pReal + endif + if (dNeq0(traction_n)) then + dudotn_dt = udotn*prm%n/traction_n + else + dudotn_dt = 0.0_pReal + endif forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + dudotd_dt*projection_d(k,l)*projection_d(m,n) & @@ -164,7 +172,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, + udott*projection_t & + udotn*projection_n enddo - + end associate end subroutine kinematics_slipplane_opening_LiAndItsTangent diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 7f7994959..658d658eb 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -11,17 +11,17 @@ module kinematics_thermal_expansion use math use lattice use material - + implicit none private - + type :: tParameters real(pReal), allocatable, dimension(:,:,:) :: & expansion end type tParameters - + type(tParameters), dimension(:), allocatable :: param - + public :: & kinematics_thermal_expansion_init, & kinematics_thermal_expansion_initialStrain, & @@ -35,36 +35,35 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine kinematics_thermal_expansion_init - - integer :: & - Ninstance, & - p, i + + integer :: Ninstance,p,i,instance real(pReal), dimension(:), allocatable :: & temp - + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'; flush(6) - + Ninstance = count(phase_kinematics == KINEMATICS_thermal_expansion_ID) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) - - do p = 1, size(phase_kinematics) + + do p = 1, size(config_phase) if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle - + associate(config => config_phase(p)) + ! ToDo: Here we need to decide how to extend the concept of instances to ! kinetics and sources. I would suggest that the same mechanism exists at maximum once per phase - + ! read up to three parameters (constant, linear, quadratic with T) - temp = config_phase(p)%getFloats('thermal_expansion11') + temp = config%getFloats('thermal_expansion11') !lattice_thermalExpansion33(1,1,1:size(temp),p) = temp - temp = config_phase(p)%getFloats('thermal_expansion22', & - defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) + temp = config%getFloats('thermal_expansion22', & + defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) !lattice_thermalExpansion33(2,2,1:size(temp),p) = temp - temp = config_phase(p)%getFloats('thermal_expansion33', & - defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) + temp = config%getFloats('thermal_expansion33', & + defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) + end associate enddo end subroutine kinematics_thermal_expansion_init @@ -74,14 +73,13 @@ end subroutine kinematics_thermal_expansion_init !> @brief report initial thermal strain based on current temperature deviation from reference !-------------------------------------------------------------------------------------------------- pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) - + integer, intent(in) :: & phase, & homog, offset real(pReal), dimension(3,3) :: & kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) - - + kinematics_thermal_expansion_initialStrain = & (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * & lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient @@ -89,15 +87,15 @@ pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * & lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient - + end function kinematics_thermal_expansion_initialStrain !-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient +!> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) - + integer, intent(in) :: & ipc, & !< grain number ip, & !< integration point number @@ -110,15 +108,15 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, phase, & homog, offset real(pReal) :: & - T, TRef, TDot - + T, TRef, TDot + phase = material_phaseAt(ipc,el) homog = material_homogenizationAt(el) offset = thermalMapping(homog)%p(ip,el) T = temperature(homog)%p(offset) TDot = temperatureRate(homog)%p(offset) TRef = lattice_referenceTemperature(phase) - + Li = TDot * ( & lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient @@ -129,8 +127,8 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & ) - dLi_dTstar = 0.0_pReal - + dLi_dTstar = 0.0_pReal + end subroutine kinematics_thermal_expansion_LiAndItsTangent end module kinematics_thermal_expansion diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index b89ec2ef7..269e3b326 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -95,23 +95,15 @@ subroutine source_damage_anisoBrittle_init 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%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)) @@ -120,8 +112,11 @@ subroutine source_damage_anisoBrittle_init 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' + 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' + 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 @@ -143,8 +138,6 @@ subroutine source_damage_anisoBrittle_init enddo - end associate - NofMyPhase=count(material_phaseAt==p) * discretization_nIP instance = source_damage_anisoBrittle_instance(p) sourceOffset = source_damage_anisoBrittle_offset(p) @@ -153,6 +146,8 @@ subroutine source_damage_anisoBrittle_init sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol source_damage_anisoBrittle_Ncleavage(1:size(param(instance)%Ncleavage),instance) = param(instance)%Ncleavage + + end associate enddo end subroutine source_damage_anisoBrittle_init diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 9e48fc343..7e8ba66cf 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -85,25 +85,20 @@ subroutine source_damage_anisoDuctile_init 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%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal) + prm%N = config%getFloat('anisoductile_ratesensitivity') + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) + prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip)) + + prm%totalNslip = sum(prm%Nslip) ! expand: family => system + prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip) - 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' !-------------------------------------------------------------------------------------------------- @@ -125,8 +120,6 @@ subroutine source_damage_anisoDuctile_init enddo - end associate - NofMyPhase=count(material_phaseAt==p) * discretization_nIP instance = source_damage_anisoDuctile_instance(p) sourceOffset = source_damage_anisoDuctile_offset(p) @@ -134,6 +127,7 @@ subroutine source_damage_anisoDuctile_init call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol + end associate enddo end subroutine source_damage_anisoDuctile_init diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 83c4286b0..dd282f920 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -80,18 +80,15 @@ subroutine source_damage_isoBrittle_init 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' @@ -115,15 +112,14 @@ subroutine source_damage_isoBrittle_init enddo + 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 + end associate - - 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 diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 59f68c5f1..48f84400f 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -77,18 +77,15 @@ subroutine source_damage_isoDuctile_init 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' @@ -112,9 +109,6 @@ subroutine source_damage_isoDuctile_init enddo - end associate - - NofMyPhase=count(material_phaseAt==p) * discretization_nIP instance = source_damage_isoDuctile_instance(p) sourceOffset = source_damage_isoDuctile_offset(p) @@ -122,6 +116,7 @@ subroutine source_damage_isoDuctile_init call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol + end associate enddo end subroutine source_damage_isoDuctile_init diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index a039dbff6..5d723f091 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -59,14 +59,16 @@ subroutine source_thermal_dissipation_init enddo if (all(phase_source(:,p) /= SOURCE_THERMAL_DISSIPATION_ID)) cycle + associate(config => config_phase(p)) instance = source_thermal_dissipation_instance(p) - param(instance)%kappa = config_phase(p)%getFloat('dissipation_coldworkcoeff') + param(instance)%kappa = config%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) + end associate enddo end subroutine source_thermal_dissipation_init @@ -75,7 +77,7 @@ end subroutine source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstances dissipation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase) +subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) integer, intent(in) :: & phase @@ -85,14 +87,14 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar Lp real(pReal), intent(out) :: & TDot, & - dTDOT_dT + 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 diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index 182ae5f43..bf1f50bf3 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -63,19 +63,21 @@ subroutine source_thermal_externalheat_init enddo if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle + associate(config => config_phase(p)) 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)%time = config%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)) + param(instance)%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(param(instance)%time)) call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) + end associate enddo end subroutine source_thermal_externalheat_init