From 431416b5ba2d57c7f29dee6b34bf1284fbfb2681 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 31 May 2021 23:02:51 +0200 Subject: [PATCH] Separating functionality lattice should become a module with static data and functions --- PRIVATE | 2 +- src/lattice.f90 | 43 ++---------- src/phase.f90 | 20 +++++- src/phase_mechanical.f90 | 81 ++++++++++------------ src/phase_mechanical_elastic.f90 | 6 +- src/phase_mechanical_plastic_dislotwin.f90 | 8 +-- src/phase_mechanical_plastic_nonlocal.f90 | 6 +- src/phase_thermal.f90 | 4 +- 8 files changed, 71 insertions(+), 99 deletions(-) diff --git a/PRIVATE b/PRIVATE index 185cb53be..5ecf7d74a 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 185cb53be76eded17565c5fa91bd9b4499cda4b8 +Subproject commit 5ecf7d74a79da155aed31e776963c7bc10eb9890 diff --git a/src/lattice.f90 b/src/lattice.f90 index d5611f3d2..d366674c7 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -365,25 +365,11 @@ module lattice 1, 1, 1, 1,-2, 1 & ],pReal),shape(BCT_SYSTEMSLIP)) !< bct slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x) - - enum, bind(c); enumerator :: & - lattice_UNDEFINED_ID, & - lattice_FCC_ID, & - lattice_BCC_ID, & - lattice_HEX_ID, & - lattice_BCT_ID - end enum - ! SHOULD NOT BE PART OF LATTICE BEGIN real(pReal), dimension(:), allocatable, public, protected :: & - lattice_mu, lattice_nu, & - lattice_mu_phi, & - lattice_rho, & - lattice_c_p + lattice_mu, lattice_nu real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_C66 - integer(kind(lattice_UNDEFINED_ID)), dimension(:), allocatable, public, protected :: & - lattice_structure ! SHOULD NOT BE PART OF LATTICE END interface lattice_forestProjection_edge @@ -396,10 +382,6 @@ module lattice public :: & lattice_init, & - lattice_FCC_ID, & - lattice_BCC_ID, & - lattice_HEX_ID, & - lattice_BCT_ID, & lattice_equivalent_nu, & lattice_equivalent_mu, & lattice_applyLatticeSymmetry33, & @@ -446,11 +428,9 @@ subroutine lattice_init phases => config_material%get('phase') Nphases = phases%length - allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID) - allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) + allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) - allocate(lattice_rho, & - lattice_mu, lattice_nu,& + allocate(lattice_mu, lattice_nu,& source=[(0.0_pReal,i=1,Nphases)]) do ph = 1, phases%length @@ -466,19 +446,6 @@ subroutine lattice_init lattice_C66(3,3,ph) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal) lattice_C66(6,6,ph) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal) - select case(phase%get_asString('lattice')) - case('cF') - lattice_structure(ph) = lattice_FCC_ID - case('cI') - lattice_structure(ph) = lattice_BCC_ID - case('hP') - lattice_structure(ph) = lattice_HEX_ID - case('tI') - lattice_structure(ph) = lattice_BCT_ID - case default - call IO_error(130,ext_msg='lattice_init: '//phase%get_asString('lattice')) - end select - lattice_C66(1:6,1:6,ph) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,ph),phase%get_asString('lattice')) lattice_nu(ph) = lattice_equivalent_nu(lattice_C66(1:6,1:6,ph),'voigt') @@ -489,8 +456,6 @@ subroutine lattice_init if (abs(lattice_C66(i,i,ph))>>'; flush(IO_STDOUT) @@ -365,9 +371,19 @@ subroutine phase_init materials => config_material%get('material') phases => config_material%get('phase') - + allocate(phase_lattice(phases%length)) + allocate(phase_cOverA(phases%length),source=-1.0_pReal) + allocate(phase_rho(phases%length)) allocate(phase_orientation0(phases%length)) + do ph = 1,phases%length + phase => phases%get(ph) + phase_lattice(ph) = phase%get_asString('lattice') + if (all(phase_lattice(ph) /= ['cF','cI','hP','tI'])) & + call IO_error(130,ext_msg='phase_init: '//phase%get_asString('lattice')) + if (any(phase_lattice(ph) == ['hP','tI'])) & + phase_cOverA(ph) = phase%get_asFloat('c/a') + phase_rho(ph) = phase%get_asFloat('rho',defaultVal=0.0_pReal) allocate(phase_orientation0(ph)%data(count(material_phaseID==ph))) enddo diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 57d1098d7..69f57fc27 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -896,55 +896,46 @@ subroutine crystallite_results(group,ph) integer :: ou real(pReal), allocatable, dimension(:,:) :: selected_rotations - character(len=:), allocatable :: structureLabel - call results_closeGroup(results_addGroup(group//'/mechanical')) + call results_closeGroup(results_addGroup(group//'/mechanical')) - do ou = 1, size(output_constituent(ph)%label) + do ou = 1, size(output_constituent(ph)%label) - select case (output_constituent(ph)%label(ou)) - case('F') - call results_writeDataset(group//'/mechanical/',phase_mechanical_F(ph)%data,'F',& - 'deformation gradient','1') - case('F_e') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Fe(ph)%data,'F_e',& - 'elastic deformation gradient','1') - case('F_p') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Fp(ph)%data,'F_p', & - 'plastic deformation gradient','1') - case('F_i') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Fi(ph)%data,'F_i', & - 'inelastic deformation gradient','1') - case('L_p') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Lp(ph)%data,'L_p', & - 'plastic velocity gradient','1/s') - case('L_i') - call results_writeDataset(group//'/mechanical/',phase_mechanical_Li(ph)%data,'L_i', & - 'inelastic velocity gradient','1/s') - case('P') - call results_writeDataset(group//'/mechanical/',phase_mechanical_P(ph)%data,'P', & - 'first Piola-Kirchhoff stress','Pa') - case('S') - call results_writeDataset(group//'/mechanical/',phase_mechanical_S(ph)%data,'S', & - 'second Piola-Kirchhoff stress','Pa') - case('O') - select case(lattice_structure(ph)) - case(lattice_FCC_ID) - structureLabel = 'cF' - case(lattice_BCC_ID) - structureLabel = 'cI' - case(lattice_HEX_ID) - structureLabel = 'hP' - case(lattice_BCT_ID) - structureLabel = 'tI' - end select - selected_rotations = select_rotations(phase_orientation(ph)%data) - call results_writeDataset(group//'/mechanical',selected_rotations,output_constituent(ph)%label(ou),& - 'crystal orientation as quaternion','q_0 (q_1 q_2 q_3)') - call results_addAttribute('lattice',structureLabel,group//'/mechanical/'//output_constituent(ph)%label(ou)) - end select - enddo + select case (output_constituent(ph)%label(ou)) + case('F') + call results_writeDataset(group//'/mechanical/',phase_mechanical_F(ph)%data,'F',& + 'deformation gradient','1') + case('F_e') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Fe(ph)%data,'F_e',& + 'elastic deformation gradient','1') + case('F_p') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Fp(ph)%data,'F_p', & + 'plastic deformation gradient','1') + case('F_i') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Fi(ph)%data,'F_i', & + 'inelastic deformation gradient','1') + case('L_p') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Lp(ph)%data,'L_p', & + 'plastic velocity gradient','1/s') + case('L_i') + call results_writeDataset(group//'/mechanical/',phase_mechanical_Li(ph)%data,'L_i', & + 'inelastic velocity gradient','1/s') + case('P') + call results_writeDataset(group//'/mechanical/',phase_mechanical_P(ph)%data,'P', & + 'first Piola-Kirchhoff stress','Pa') + case('S') + call results_writeDataset(group//'/mechanical/',phase_mechanical_S(ph)%data,'S', & + 'second Piola-Kirchhoff stress','Pa') + case('O') + selected_rotations = select_rotations(phase_orientation(ph)%data) + call results_writeDataset(group//'/mechanical',selected_rotations,output_constituent(ph)%label(ou),& + 'crystal orientation as quaternion','q_0 (q_1 q_2 q_3)') + call results_addAttribute('lattice',phase_lattice(ph),group//'/mechanical/'//output_constituent(ph)%label(ou)) + if (any(phase_lattice(ph) == ['hP', 'tI'])) & + call results_addAttribute('c/a',phase_cOverA(ph),group//'/mechanical/'//output_constituent(ph)%label(ou)) + end select + enddo contains diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 94622d359..53ed22613 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -4,7 +4,7 @@ submodule(phase:mechanical) elastic real(pReal), dimension(6,6) :: & C66 !< Elastic constants in Voig notation end type tParameters - + type(tParameters), allocatable, dimension(:) :: param contains @@ -30,13 +30,13 @@ module subroutine elastic_init(phases) print'(a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT) allocate(param(phases%length)) - + do ph = 1, phases%length phase => phases%get(ph) mech => phase%get('mechanical') elastic => mech%get('elastic') if (elastic%get_asString('type') /= 'Hooke') call IO_error(200,ext_msg=elastic%get_asString('type')) - + associate(prm => param(ph)) struct = phase%get_asString('lattice') if (struct /= 'cI' .and. struct /= 'cF' .and. struct /= 'hP' .and. struct /= 'tI') & diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index fdfd630c6..254dda56e 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -202,7 +202,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%n0_sl = lattice_slip_normal(N_sl,phase%get_asString('lattice'),& phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - prm%fccTwinTransNucleation = lattice_structure(ph) == lattice_FCC_ID .and. (N_sl(1) == 12) + prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. (N_sl(1) == 12) if(prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl)) @@ -231,7 +231,7 @@ module function plastic_dislotwin_init() result(myPlasticity) ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) ! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 prm%omega = pl%get_asFloat('omega', defaultVal = 1000.0_pReal) & - * merge(12.0_pReal,8.0_pReal,any(lattice_structure(ph) == [lattice_FCC_ID,lattice_HEX_ID])) + * merge(12.0_pReal,8.0_pReal,any(phase_lattice(ph) == ['cF','hP'])) ! expand: family => system rho_mob_0 = math_expand(rho_mob_0, N_sl) @@ -340,7 +340,7 @@ module function plastic_dislotwin_init() result(myPlasticity) pl%get_asFloat('a_cI', defaultVal=0.0_pReal), & pl%get_asFloat('a_cF', defaultVal=0.0_pReal)) - if (lattice_structure(ph) /= lattice_FCC_ID) then + if (phase_lattice(ph) /= 'cF') then prm%dot_N_0_tr = pl%get_as1dFloat('dot_N_0_tr') prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr) endif @@ -355,7 +355,7 @@ module function plastic_dislotwin_init() result(myPlasticity) if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr' if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr' if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr' - if (lattice_structure(ph) /= lattice_FCC_ID) then + if (phase_lattice(ph) /= 'cF') then if (any(prm%dot_N_0_tr < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tr' endif else transActive diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 960e668fa..7ece6f4ed 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -622,7 +622,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) ! coefficients are corrected for the line tension effect ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) - if (any(lattice_structure(ph) == [LATTICE_bcc_ID,LATTICE_fcc_ID])) then + if (any(phase_lattice(ph) == ['cI','cF'])) then myInteractionMatrix = prm%h_sl_sl & * spread(( 1.0_pReal - prm%f_F & + prm%f_F & @@ -1055,7 +1055,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & !**************************************************************************** !*** dislocation multiplication rhoDotMultiplication = 0.0_pReal - isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) then + isBCC: if (phase_lattice(ph) == 'cI') then forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path @@ -1111,7 +1111,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & + rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent ! annihilated screw dipoles leave edge jogs behind on the colinear system - if (lattice_structure(ph) == LATTICE_fcc_ID) & + if (phase_lattice(ph) == 'cF') & forall (s = 1:ns, prm%colinearSystem(s) > 0) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & * 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 141cf0249..808cd7cf2 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -196,8 +196,8 @@ module function phase_mu_T(co,ce) result(mu) real(pReal) :: mu - mu = lattice_rho(material_phaseID(co,ce)) & - * param(material_phaseID(co,ce))%C_p + mu = phase_rho(material_phaseID(co,ce)) & + * param(material_phaseID(co,ce))%C_p end function phase_mu_T