From 4d227fab2b3c789b068b6432e58a417f2926d53e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Feb 2020 14:34:19 +0100 Subject: [PATCH] polishing --- src/kinematics_cleavage_opening.f90 | 1 + src/kinematics_slipplane_opening.f90 | 7 +- src/kinematics_thermal_expansion.f90 | 4 +- src/lattice.f90 | 120 ++++++++++++++------------- src/source_damage_anisoBrittle.f90 | 11 ++- src/source_damage_isoBrittle.f90 | 12 ++- src/source_damage_isoDuctile.f90 | 10 +-- src/source_thermal_dissipation.f90 | 1 + src/source_thermal_externalheat.f90 | 8 +- 9 files changed, 88 insertions(+), 86 deletions(-) diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index a68606f07..1a8095ff4 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -110,6 +110,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + integer :: & homog, damageOffset, & i, k, l, m, n diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index ad1fd2c20..a25d110e4 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -29,9 +29,9 @@ module kinematics_slipplane_opening real(pReal), dimension(:), allocatable :: & critLoad real(pReal), dimension(:,:), allocatable :: & - slip_direction, & - slip_normal, & - slip_transverse + slip_direction, & + slip_normal, & + slip_transverse end type tParameters type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) @@ -113,6 +113,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + real(pReal), dimension(3,3) :: & projection_d, projection_t, projection_n !< projection modes 3x3 tensor integer :: & diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 822d9ca02..72063e56e 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -41,8 +41,7 @@ contains subroutine kinematics_thermal_expansion_init integer :: Ninstance,p,i - real(pReal), dimension(:), allocatable :: & - temp + real(pReal), dimension(:), allocatable :: temp write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'; flush(6) @@ -88,6 +87,7 @@ pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) phase, & homog, & offset + real(pReal), dimension(3,3) :: & kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) diff --git a/src/lattice.f90 b/src/lattice.f90 index bc074056b..11576d7ca 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -463,7 +463,7 @@ subroutine lattice_init allocate(lattice_thermalConductivity (3,3,Nphases), source=0.0_pReal) allocate(lattice_damageDiffusion (3,3,Nphases), source=0.0_pReal) - + allocate(lattice_damageMobility,& lattice_massDensity,lattice_specificHeat, & lattice_mu, lattice_nu,& @@ -500,7 +500,7 @@ subroutine lattice_init call IO_error(130,ext_msg='lattice_init: '//trim(structure)) end select - lattice_C66(1:6,1:6,p) = symmetrizeC66(lattice_structure(p),lattice_C66(1:6,1:6,p)) + lattice_C66(1:6,1:6,p) = symmetrizeC66(lattice_C66(1:6,1:6,p),structure) lattice_mu(p) = 0.2_pReal *(lattice_C66(1,1,p) -lattice_C66(1,2,p) +3.0_pReal*lattice_C66(4,4,p)) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 lattice_nu(p) = ( lattice_C66(1,1,p) +4.0_pReal*lattice_C66(1,2,p) -2.0_pReal*lattice_C66(4,4,p)) & @@ -527,9 +527,9 @@ subroutine lattice_init lattice_DamageDiffusion(2,2,p) = config_phase(p)%getFloat('damage_diffusion22',defaultVal=0.0_pReal) lattice_DamageDiffusion(3,3,p) = config_phase(p)%getFloat('damage_diffusion33',defaultVal=0.0_pReal) lattice_DamageDiffusion(1:3,1:3,p) = lattice_symmetrize33(lattice_DamageDiffusion(1:3,1:3,p),structure) - + lattice_DamageMobility(p) = config_phase(p)%getFloat( 'damage_mobility',defaultVal=0.0_pReal) - + enddo end subroutine lattice_init @@ -690,7 +690,7 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, & C_target_unrotated66(1,3) = C_bar66(1,3) C_target_unrotated66(3,3) = C_bar66(3,3) C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2))) - C_target_unrotated66 = symmetrizeC66(LATTICE_HEX_ID,C_target_unrotated66) + C_target_unrotated66 = symmetrizeC66(C_target_unrotated66,'hex') elseif (structure_target(1:3) == 'bcc') then if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) & call IO_error(134,ext_msg='lattice_C66_trans: '//trim(structure_target)) @@ -1797,7 +1797,7 @@ function lattice_symmetrize33(T,structure) result(T_sym) T_sym(2,2) = T(2,2) T_sym(3,3) = T(3,3) case default - call IO_error(137,ext_msg='lattice_symmetrize33: '//trim(structure)) + call IO_error(137,ext_msg='lattice_symmetrize33: '//trim(structure)) end select end function lattice_symmetrize33 @@ -1807,73 +1807,75 @@ end function lattice_symmetrize33 !> @brief Symmetrizes stiffness matrix according to lattice type !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 !-------------------------------------------------------------------------------------------------- -pure function symmetrizeC66(struct,C66) +function symmetrizeC66(C66,structure) result(C66_sym) + + real(pReal), dimension(6,6) :: C66_sym - integer(kind(LATTICE_undefined_ID)), intent(in) :: struct real(pReal), dimension(6,6), intent(in) :: C66 - real(pReal), dimension(6,6) :: symmetrizeC66 + character(len=*), intent(in) :: structure + integer :: j,k - symmetrizeC66 = 0.0_pReal + C66_sym = 0.0_pReal - select case(struct) - case (LATTICE_iso_ID) + select case(structure) + case ('iso') do k=1,3 do j=1,3 - symmetrizeC66(k,j) = C66(1,2) + C66_sym(k,j) = C66(1,2) enddo - symmetrizeC66(k,k) = C66(1,1) - symmetrizeC66(k+3,k+3) = 0.5_pReal*(C66(1,1)-C66(1,2)) + C66_sym(k,k) = C66(1,1) + C66_sym(k+3,k+3) = 0.5_pReal*(C66(1,1)-C66(1,2)) enddo - case (LATTICE_fcc_ID,LATTICE_bcc_ID) + case ('fcc','bcc') do k=1,3 do j=1,3 - symmetrizeC66(k,j) = C66(1,2) + C66_sym(k,j) = C66(1,2) enddo - symmetrizeC66(k,k) = C66(1,1) - symmetrizeC66(k+3,k+3) = C66(4,4) + C66_sym(k,k) = C66(1,1) + C66_sym(k+3,k+3) = C66(4,4) enddo - case (LATTICE_hex_ID) - symmetrizeC66(1,1) = C66(1,1) - symmetrizeC66(2,2) = C66(1,1) - symmetrizeC66(3,3) = C66(3,3) - symmetrizeC66(1,2) = C66(1,2) - symmetrizeC66(2,1) = C66(1,2) - symmetrizeC66(1,3) = C66(1,3) - symmetrizeC66(3,1) = C66(1,3) - symmetrizeC66(2,3) = C66(1,3) - symmetrizeC66(3,2) = C66(1,3) - symmetrizeC66(4,4) = C66(4,4) - symmetrizeC66(5,5) = C66(4,4) - symmetrizeC66(6,6) = 0.5_pReal*(C66(1,1)-C66(1,2)) - case (LATTICE_ort_ID) - symmetrizeC66(1,1) = C66(1,1) - symmetrizeC66(2,2) = C66(2,2) - symmetrizeC66(3,3) = C66(3,3) - symmetrizeC66(1,2) = C66(1,2) - symmetrizeC66(2,1) = C66(1,2) - symmetrizeC66(1,3) = C66(1,3) - symmetrizeC66(3,1) = C66(1,3) - symmetrizeC66(2,3) = C66(2,3) - symmetrizeC66(3,2) = C66(2,3) - symmetrizeC66(4,4) = C66(4,4) - symmetrizeC66(5,5) = C66(5,5) - symmetrizeC66(6,6) = C66(6,6) - case (LATTICE_bct_ID) - symmetrizeC66(1,1) = C66(1,1) - symmetrizeC66(2,2) = C66(1,1) - symmetrizeC66(3,3) = C66(3,3) - symmetrizeC66(1,2) = C66(1,2) - symmetrizeC66(2,1) = C66(1,2) - symmetrizeC66(1,3) = C66(1,3) - symmetrizeC66(3,1) = C66(1,3) - symmetrizeC66(2,3) = C66(1,3) - symmetrizeC66(3,2) = C66(1,3) - symmetrizeC66(4,4) = C66(4,4) - symmetrizeC66(5,5) = C66(4,4) - symmetrizeC66(6,6) = C66(6,6) + case ('hex') + C66_sym(1,1) = C66(1,1) + C66_sym(2,2) = C66(1,1) + C66_sym(3,3) = C66(3,3) + C66_sym(1,2) = C66(1,2) + C66_sym(2,1) = C66(1,2) + C66_sym(1,3) = C66(1,3) + C66_sym(3,1) = C66(1,3) + C66_sym(2,3) = C66(1,3) + C66_sym(3,2) = C66(1,3) + C66_sym(4,4) = C66(4,4) + C66_sym(5,5) = C66(4,4) + C66_sym(6,6) = 0.5_pReal*(C66(1,1)-C66(1,2)) + case ('ort') + C66_sym(1,1) = C66(1,1) + C66_sym(2,2) = C66(2,2) + C66_sym(3,3) = C66(3,3) + C66_sym(1,2) = C66(1,2) + C66_sym(2,1) = C66(1,2) + C66_sym(1,3) = C66(1,3) + C66_sym(3,1) = C66(1,3) + C66_sym(2,3) = C66(2,3) + C66_sym(3,2) = C66(2,3) + C66_sym(4,4) = C66(4,4) + C66_sym(5,5) = C66(5,5) + C66_sym(6,6) = C66(6,6) + case ('bct') + C66_sym(1,1) = C66(1,1) + C66_sym(2,2) = C66(1,1) + C66_sym(3,3) = C66(3,3) + C66_sym(1,2) = C66(1,2) + C66_sym(2,1) = C66(1,2) + C66_sym(1,3) = C66(1,3) + C66_sym(3,1) = C66(1,3) + C66_sym(2,3) = C66(1,3) + C66_sym(3,2) = C66(1,3) + C66_sym(4,4) = C66(4,4) + C66_sym(5,5) = C66(4,4) + C66_sym(6,6) = C66(6,6) case default - symmetrizeC66 = C66 + call IO_error(137,ext_msg='symmetrizeC66: '//trim(structure)) end select end function symmetrizeC66 diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 4e394a24b..91c90adc2 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -22,7 +22,7 @@ module source_damage_anisoBrittle source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? source_damage_anisoBrittle_instance !< instance of source mechanism - type :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & aTol, & sdot_0, & @@ -166,11 +166,10 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & - + prm%sdot_0* & - ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%N + & - (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%N + & - (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%N)/ & - prm%critDisp(i) + + prm%sdot_0 / prm%critDisp(i) & + * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%N + & + (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%N + & + (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%N) enddo end associate diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index a31202273..8856dc4e0 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -125,9 +125,8 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) real(pReal) :: & strainenergy - phase = material_phaseAt(ipc,el) !< phase ID at ipc,ip,el - constituent = material_phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el - ! ToDo: capability for multiple instances of SAME source within given phase. Needs Ninstance loop from here on! + phase = material_phaseAt(ipc,el) !< phase ID at ipc,ip,el + constituent = material_phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el sourceOffset = source_damage_isoBrittle_offset(phase) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) @@ -169,10 +168,9 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD sourceOffset = source_damage_isoBrittle_offset(phase) associate(prm => param(source_damage_isoBrittle_instance(phase))) - localphiDot = (1.0_pReal - phi)**(prm%n - 1.0_pReal) - & - phi*sourceState(phase)%p(sourceOffset)%state(1,constituent) - dLocalphiDot_dPhi = - (prm%n - 1.0_pReal)* & - (1.0_pReal - phi)**max(0.0_pReal,prm%n - 2.0_pReal) & + localphiDot = (1.0_pReal - phi)**(prm%n - 1.0_pReal) & + - phi*sourceState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = - (prm%n - 1.0_pReal)* (1.0_pReal - phi)**max(0.0_pReal,prm%n - 2.0_pReal) & - sourceState(phase)%p(sourceOffset)%state(1,constituent) end associate diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index edfe8917f..2189d5b2b 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -17,10 +17,10 @@ module source_damage_isoDuctile private integer, dimension(:), allocatable :: & - source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? - source_damage_isoDuctile_instance !< instance of damage source mechanism + source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? + source_damage_isoDuctile_instance !< instance of damage source mechanism - type, private :: tParameters !< container type for internal constitutive parameters + type, private :: tParameters !< container type for internal constitutive parameters real(pReal) :: & critPlasticStrain, & N, & @@ -29,7 +29,7 @@ module source_damage_isoDuctile output end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) public :: & @@ -126,7 +126,7 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el) associate(prm => param(source_damage_isoDuctile_instance(phase))) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain + sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain ! ToDo: abs for slip rate? end associate end subroutine source_damage_isoDuctile_dotState diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index acaf30068..c03d66c24 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -86,6 +86,7 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar Tstar real(pReal), intent(in), dimension(3,3) :: & Lp + real(pReal), intent(out) :: & TDot, & dTDot_dT diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index 0380a2be3..6cb8cdf6e 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -122,15 +122,15 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas sourceOffset = source_thermal_externalheat_offset(phase) associate(prm => param(source_thermal_externalheat_instance(phase))) - do interval = 1, prm%nIntervals ! scan through all rate segments + do interval = 1, prm%nIntervals ! scan through all rate segments frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%time(interval)) & - / (prm%time(interval+1) - prm%time(interval)) ! fractional time within segment + / (prm%time(interval+1) - prm%time(interval)) ! fractional time within segment if ( (frac_time < 0.0_pReal .and. interval == 1) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & TDot = prm%heat_rate(interval ) * (1.0_pReal - frac_time) + & - prm%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries... - ! ...or extrapolate if outside of bounds + prm%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries... + ! ...or extrapolate if outside of bounds enddo dTDot_dT = 0.0 end associate