From 37d511062d70fdfedf89444a0e68270ae6adc313 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 09:24:00 +0200 Subject: [PATCH 01/18] simplified (using new mapping) --- src/phase_mechanical.f90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 29e9e8ada..e2ad8f968 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -971,11 +971,9 @@ subroutine crystallite_results(group,ph) integer, intent(in) :: ph type(rotation), dimension(:,:,:), intent(in) :: dataset - real(pReal), allocatable, dimension(:,:) :: select_rotations + real(pReal), dimension(4,count(material_phaseID==ph)) :: select_rotations integer :: el,ip,co,j - allocate(select_rotations(4,count(material_phaseAt==ph)*homogenization_maxNconstituents*discretization_nIPs)) - j=0 do el = 1, size(material_phaseAt,2) do ip = 1, discretization_nIPs From 00fcf069fbd94756907e2d3dd8255fafabf52a92 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 09:33:58 +0200 Subject: [PATCH 02/18] clearer names and reporting --- src/material.f90 | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index add9fd1c3..e072d9a1c 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -18,14 +18,13 @@ module material integer, dimension(:), allocatable, public, protected :: & homogenization_Nconstituents !< number of grains in each homogenization + integer, public, protected :: & + homogenization_maxNconstituents !< max number of grains in any homogenization character(len=:), public, protected, allocatable, dimension(:) :: & material_name_phase, & !< name of each phase material_name_homogenization !< name of each homogenization - integer, public, protected :: & - homogenization_maxNconstituents !< max number of grains in any USED homogenization - integer, dimension(:), allocatable, public, protected :: & ! (elem) material_homogenizationAt, & !< homogenization ID of each element material_homogenizationID, & !< per cell @@ -55,8 +54,8 @@ subroutine material_init(restart) print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) - call material_parseMaterial - print*, 'Material parsed' + call parse + print*, 'parsed material.yaml' if (.not. restart) then @@ -69,11 +68,10 @@ subroutine material_init(restart) end subroutine material_init - !-------------------------------------------------------------------------------------------------- -!> @brief parses the material part in the material configuration file +!> @brief Parse material.yaml to get the global structure !-------------------------------------------------------------------------------------------------- -subroutine material_parseMaterial +subroutine parse() class(tNode), pointer :: materials, & !> list of materials material, & !> material definition @@ -155,7 +153,7 @@ subroutine material_parseMaterial enddo -end subroutine material_parseMaterial +end subroutine parse !-------------------------------------------------------------------------------------------------- From 77d1ed465e92db3a465ebe3640f94ec0c216a47f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 10:20:00 +0200 Subject: [PATCH 03/18] divide and test --- src/IO.f90 | 57 ++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 43 insertions(+), 14 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 81745441a..8f94a40c7 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -20,6 +20,9 @@ module IO character, parameter, public :: & IO_EOL = new_line('DAMASK'), & !< end of line character IO_COMMENT = '#' + character, parameter :: & + CR = achar(13), & + LF = IO_EOL character(len=*), parameter :: & IO_DIVIDER = '───────────────────'//& '───────────────────'//& @@ -112,8 +115,8 @@ end function IO_readlines !-------------------------------------------------------------------------------------------------- -!> @brief Read whole file. -!> @details ensures that the string ends with a new line (expected UNIX behavior) +!> @brief Read ASCII file. +!> @details Proper Unix style (LF line endings and LF at EOF) is ensured. !-------------------------------------------------------------------------------------------------- function IO_read(fileName) result(fileContent) @@ -124,7 +127,6 @@ function IO_read(fileName) result(fileContent) fileLength, & fileUnit, & myStat - character, parameter :: CR = achar(13) inquire(file = fileName, size=fileLength) @@ -141,16 +143,8 @@ function IO_read(fileName) result(fileContent) if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName)) close(fileUnit) - foundCRLF: if (scan(fileContent(:index(fileContent,IO_EOL)),CR) /= 0) then - CRLF2LF: block - integer :: c - do c=1, len(fileContent) - if (fileContent(c:c) == CR) fileContent(c:c) = ' ' - enddo - end block CRLF2LF - endif foundCRLF - - if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF + if (scan(fileContent(:index(fileContent,LF)),CR//LF) /= 0) fileContent = CRLF2LF(fileContent) + if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF end function IO_read @@ -633,10 +627,36 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) end subroutine IO_warning +!-------------------------------------------------------------------------------------------------- +!> @brief Convert Windows (CRLF) to Unix (LF) line endings. +!-------------------------------------------------------------------------------------------------- +pure function CRLF2LF(string) + + character(len=*), intent(in) :: string + character(len=:), allocatable :: CRLF2LF + + integer :: c,n + + + allocate(character(len=len_trim(string))::CRLF2LF) + if (len(CRLF2LF) == 0) return + + n = 0 + do c=1, len_trim(string) + CRLF2LF(c-n:c-n) = string(c:c) + if (c == len_trim(string)) exit + if (string(c:c+1) == CR//LF) n = n + 1 + enddo + + CRLF2LF = CRLF2LF(:c-n) + +end function + + !-------------------------------------------------------------------------------------------------- !> @brief Check correctness of some IO functions. !-------------------------------------------------------------------------------------------------- -subroutine selfTest +subroutine selfTest() integer, dimension(:), allocatable :: chunkPos character(len=:), allocatable :: str @@ -671,6 +691,15 @@ subroutine selfTest chunkPos = IO_stringPos(str) if(3112019 /= IO_intValue(str,chunkPos,2)) error stop 'IO_intValue' + if (CRLF2LF('') /= '') error stop 'CRLF2LF/0' + if (CRLF2LF(LF) /= LF) error stop 'CRLF2LF/1a' + if (CRLF2LF(CR//LF) /= LF) error stop 'CRLF2LF/1b' + if (CRLF2LF(' '//LF) /= ' '//LF) error stop 'CRLF2LF/2a' + if (CRLF2LF(' '//CR//LF) /= ' '//LF) error stop 'CRLF2LF/2b' + if (CRLF2LF('A'//CR//LF//'B') /= 'A'//LF//'B') error stop 'CRLF2LF/3' + if (CRLF2LF('A'//CR//LF//'B'//CR//LF) /= & + 'A'//LF//'B'//LF) error stop 'CRLF2LF/4' + if(.not. IO_isBlank(' ')) error stop 'IO_isBlank/1' if(.not. IO_isBlank(' #isBlank')) error stop 'IO_isBlank/2' if( IO_isBlank(' i#s')) error stop 'IO_isBlank/3' From d2855913b55ffbd4f7c2e71732513dd9ba79c7af Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 11:26:19 +0200 Subject: [PATCH 04/18] new mappings --- src/discretization.f90 | 2 +- src/grid/DAMASK_grid.f90 | 2 +- src/phase.f90 | 7 +++---- src/phase_damage.f90 | 11 +++++----- src/phase_mechanical.f90 | 36 +++++++++++++++----------------- src/phase_mechanical_plastic.f90 | 6 +++--- 6 files changed, 30 insertions(+), 34 deletions(-) diff --git a/src/discretization.f90 b/src/discretization.f90 index 2b19cb20b..22abbd644 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -15,7 +15,7 @@ module discretization discretization_Nelems integer, public, protected, dimension(:), allocatable :: & - discretization_materialAt + discretization_materialAt !ToDo: discretization_materialID real(pReal), public, protected, dimension(:,:), allocatable :: & discretization_IPcoords0, & diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 991c39024..778cc855e 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -284,7 +284,7 @@ program DAMASK_grid if (loadCases(l)%f_restart < 1) errorID = 839 if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then - print'(a)', ' r: 1 (constant step widths)' + print'(a)', ' r: 1 (constant step width)' else print'(a,f0.3)', ' r: ', loadCases(l)%r endif diff --git a/src/phase.f90 b/src/phase.f90 index 927ba6b2b..f0b5c4e22 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -220,12 +220,11 @@ module phase end function thermal_stress - module function integrateDamageState(dt,co,ip,el) result(broken) + module function integrateDamageState(dt,co,ce) result(broken) real(pReal), intent(in) :: dt integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + ce, & + co logical :: broken end function integrateDamageState diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 5a1d5970b..5b9d8fd5f 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -171,13 +171,12 @@ end function phase_f_phi !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -module function integrateDamageState(dt,co,ip,el) result(broken) +module function integrateDamageState(dt,co,ce) result(broken) real(pReal), intent(in) :: dt integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + ce, & + co logical :: broken integer :: & @@ -193,8 +192,8 @@ module function integrateDamageState(dt,co,ip,el) result(broken) logical :: & converged_ - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + ph = material_phaseID(co,ce) + me = material_phaseEntry(co,ce) if (damageState(ph)%sizeState == 0) then broken = .false. diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index e2ad8f968..5c547df53 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -192,8 +192,6 @@ module subroutine mechanical_init(materials,phases) phases integer :: & - el, & - ip, & co, & ce, & ph, & @@ -257,14 +255,14 @@ module subroutine mechanical_init(materials,phases) #endif enddo - do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - material => materials%get(discretization_materialAt(el)) + do ce = 1, size(discretization_materialAt,1) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) + material => materials%get(discretization_materialAt(ce)) constituents => material%get('constituents') constituent => constituents%get(co) - ph = material_phaseAt(co,el) - en = material_phaseMemberAt(co,ip,el) + ph = material_phaseID(co,ce) + en = material_phaseEntry(co,ce) call material_orientation0(co,ph,en)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) @@ -281,7 +279,7 @@ module subroutine mechanical_init(materials,phases) phase_mechanical_F(ph)%data(1:3,1:3,en) = phase_mechanical_F0(ph)%data(1:3,1:3,en) enddo - enddo; enddo + enddo ! initialize elasticity @@ -427,8 +425,8 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) broken = .true. - ph = material_phaseAt(co,el) - en = material_phaseMemberAt(co,ip,el) + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) call plastic_dependentState(co,ip,el) @@ -603,8 +601,8 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul dotState - ph = material_phaseAt(co,el) - en = material_phaseMemberAt(co,ip,el) + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) return @@ -688,8 +686,8 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res sizeDotState - ph = material_phaseAt(co,el) - en = material_phaseMemberAt(co,ip,el) + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) return @@ -728,8 +726,8 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic - ph = material_phaseAt(co,el) - en = material_phaseMemberAt(co,ip,el) + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) return @@ -846,8 +844,8 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - ph = material_phaseAt(co,el) - en = material_phaseMemberAt(co,ip,el) + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) broken = plastic_dotState(Delta_t,co,ip,el,ph,en) if(broken) return @@ -1105,7 +1103,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subF = subF0 & + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) - converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) + converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,(el-1)*discretization_nIPs + ip) endif enddo cutbackLooping diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 33dfd681f..6a053e12d 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -360,10 +360,10 @@ module subroutine plastic_dependentState(co, ip, el) en - ph = material_phaseAt(co,el) - en = material_phasememberAt(co,ip,el) + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) - plasticType: select case (phase_plasticity(material_phaseAt(co,el))) + plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_DISLOTWIN_ID) plasticType call dislotwin_dependentState(thermal_T(ph,en),ph,en) From fc4bdfb3749bacab0d7184c9be49ca8a4eb4943f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 12:03:13 +0200 Subject: [PATCH 05/18] simplified --- src/phase_mechanical_plastic_nonlocal.f90 | 82 ++++++++++------------- 1 file changed, 37 insertions(+), 45 deletions(-) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 05ba371ce..8eb90f547 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -623,7 +623,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(material_phaseAt(1,el)) == [LATTICE_bcc_ID,LATTICE_fcc_ID])) then + if (any(lattice_structure(ph) == [LATTICE_bcc_ID,LATTICE_fcc_ID])) then myInteractionMatrix = prm%h_sl_sl & * spread(( 1.0_pReal - prm%f_F & + prm%f_F & @@ -644,7 +644,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) !################################################################################################# rho0 = getRho0(ph,en) - if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then + if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,en)) @@ -1238,7 +1238,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) rhoDotFlux = 0.0_pReal - if (.not. phase_localPlasticity(material_phaseAt(1,el))) then + if (.not. phase_localPlasticity(ph)) then !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... @@ -1579,21 +1579,20 @@ end subroutine plastic_nonlocal_results !-------------------------------------------------------------------------------------------------- !> @brief populates the initial dislocation density !-------------------------------------------------------------------------------------------------- -subroutine stateInit(ini,phase,Nmembers) +subroutine stateInit(ini,phase,Nentries) type(tInitialParameters) :: & ini integer,intent(in) :: & phase, & - Nmembers + Nentries + integer :: & - i, & e, & f, & from, & upto, & - s, & - phasemember + s real(pReal), dimension(2) :: & noise, & rnd @@ -1602,49 +1601,42 @@ subroutine stateInit(ini,phase,Nmembers) totalVolume, & densityBinning, & minimumIpVolume - real(pReal), dimension(Nmembers) :: & - volume associate(stt => state(phase)) - if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume - do e = 1,discretization_Nelems - do i = 1,discretization_nIPs - if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e) - enddo - enddo - totalVolume = sum(volume) - minimumIPVolume = minval(volume) - densityBinning = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pReal / 3.0_pReal) + if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume + totalVolume = sum(geom(phase)%V_0) + minimumIPVolume = minval(geom(phase)%V_0) + densityBinning = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pReal / 3.0_pReal) - ! fill random material points with dislocation segments until the desired overall density is reached - meanDensity = 0.0_pReal - do while(meanDensity < ini%random_rho_u) - call random_number(rnd) - phasemember = nint(rnd(1)*real(Nmembers,pReal) + 0.5_pReal) - s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) - meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume - stt%rhoSglMobile(s,phasemember) = densityBinning - enddo - else ! homogeneous distribution with noise - do e = 1, Nmembers - do f = 1,size(ini%N_sl,1) - from = 1 + sum(ini%N_sl(1:f-1)) - upto = sum(ini%N_sl(1:f)) - do s = from,upto - noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), & - math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)] - stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1) - stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) - stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2) - stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2) - enddo - stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f) - stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f) + ! fill random material points with dislocation segments until the desired overall density is reached + meanDensity = 0.0_pReal + do while(meanDensity < ini%random_rho_u) + call random_number(rnd) + e = nint(rnd(1)*real(Nentries,pReal) + 0.5_pReal) + s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) + meanDensity = meanDensity + densityBinning * geom(phase)%V_0(e) / totalVolume + stt%rhoSglMobile(s,e) = densityBinning enddo - enddo - endif + else ! homogeneous distribution with noise + do e = 1, Nentries + do f = 1,size(ini%N_sl,1) + from = 1 + sum(ini%N_sl(1:f-1)) + upto = sum(ini%N_sl(1:f)) + do s = from,upto + noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), & + math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)] + stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1) + stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) + stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2) + stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2) + enddo + stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f) + stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f) + enddo + enddo + endif end associate From 212c4296cf3f7d5cfdc0009201a7e42b86cc77e6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 16:42:06 +0200 Subject: [PATCH 06/18] need to loop over ip, el --- src/phase_mechanical.f90 | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 5c547df53..cbaf149ec 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -192,6 +192,8 @@ module subroutine mechanical_init(materials,phases) phases integer :: & + el, & + ip, & co, & ce, & ph, & @@ -255,14 +257,14 @@ module subroutine mechanical_init(materials,phases) #endif enddo - do ce = 1, size(discretization_materialAt,1) - do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) - material => materials%get(discretization_materialAt(ce)) + do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + material => materials%get(discretization_materialAt(el)) constituents => material%get('constituents') constituent => constituents%get(co) - ph = material_phaseID(co,ce) - en = material_phaseEntry(co,ce) + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) call material_orientation0(co,ph,en)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) @@ -279,7 +281,7 @@ module subroutine mechanical_init(materials,phases) phase_mechanical_F(ph)%data(1:3,1:3,en) = phase_mechanical_F0(ph)%data(1:3,1:3,en) enddo - enddo + enddo; enddo ! initialize elasticity From 6d31f77debdb14de9bc2e8d4c7768d9f59a9a32e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 17:21:07 +0200 Subject: [PATCH 07/18] natural data structure --- src/discretization.f90 | 30 +++++++++++------------ src/material.f90 | 21 +++++++++++++++- src/phase.f90 | 28 ++++++++++++++++++--- src/phase_mechanical.f90 | 11 ++------- src/phase_mechanical_plastic_nonlocal.f90 | 4 +-- 5 files changed, 63 insertions(+), 31 deletions(-) diff --git a/src/discretization.f90 b/src/discretization.f90 index 22abbd644..cb74056cc 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -9,20 +9,20 @@ module discretization implicit none private - + integer, public, protected :: & discretization_nIPs, & discretization_Nelems - - integer, public, protected, dimension(:), allocatable :: & - discretization_materialAt !ToDo: discretization_materialID - real(pReal), public, protected, dimension(:,:), allocatable :: & + integer, public, protected, dimension(:), allocatable :: & + discretization_materialAt !ToDo: discretization_ID_material + + real(pReal), public, protected, dimension(:,:), allocatable :: & discretization_IPcoords0, & discretization_IPcoords, & discretization_NodeCoords0, & discretization_NodeCoords - + integer :: & discretization_sharedNodesBegin @@ -33,7 +33,7 @@ module discretization discretization_setNodeCoords contains - + !-------------------------------------------------------------------------------------------------- !> @brief stores the relevant information in globally accesible variables !-------------------------------------------------------------------------------------------------- @@ -54,20 +54,20 @@ subroutine discretization_init(materialAt,& discretization_Nelems = size(materialAt,1) discretization_nIPs = size(IPcoords0,2)/discretization_Nelems - discretization_materialAt = materialAt + discretization_materialAt = materialAt discretization_IPcoords0 = IPcoords0 discretization_IPcoords = IPcoords0 discretization_NodeCoords0 = NodeCoords0 discretization_NodeCoords = NodeCoords0 - + if(present(sharedNodesBegin)) then discretization_sharedNodesBegin = sharedNodesBegin else discretization_sharedNodesBegin = size(discretization_NodeCoords0,2) endif - + end subroutine discretization_init @@ -77,13 +77,13 @@ end subroutine discretization_init subroutine discretization_results real(pReal), dimension(:,:), allocatable :: u - + call results_closeGroup(results_addGroup('current/geometry')) - + u = discretization_NodeCoords (1:3,:discretization_sharedNodesBegin) & - discretization_NodeCoords0(1:3,:discretization_sharedNodesBegin) call results_writeDataset('current/geometry',u,'u_n','displacements of the nodes','m') - + u = discretization_IPcoords & - discretization_IPcoords0 call results_writeDataset('current/geometry',u,'u_p','displacements of the materialpoints (cell centers)','m') @@ -97,7 +97,7 @@ end subroutine discretization_results subroutine discretization_setIPcoords(IPcoords) real(pReal), dimension(:,:), intent(in) :: IPcoords - + discretization_IPcoords = IPcoords end subroutine discretization_setIPcoords @@ -109,7 +109,7 @@ end subroutine discretization_setIPcoords subroutine discretization_setNodeCoords(NodeCoords) real(pReal), dimension(:,:), intent(in) :: NodeCoords - + discretization_NodeCoords = NodeCoords end subroutine discretization_setNodeCoords diff --git a/src/material.f90 b/src/material.f90 index e072d9a1c..c62442fcd 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -16,6 +16,12 @@ module material implicit none private + type :: tRotationContainer + type(Rotation), dimension(:), allocatable :: data + end type + + type(tRotationContainer), dimension(:), allocatable :: material_orientation0_2 + integer, dimension(:), allocatable, public, protected :: & homogenization_Nconstituents !< number of grains in each homogenization integer, public, protected :: & @@ -39,6 +45,8 @@ module material material_phaseMemberAt !< position of the element within its phase instance public :: & + tRotationContainer, & + material_orientation0_2, & material_init contains @@ -88,7 +96,7 @@ subroutine parse() real(pReal) :: & frac integer :: & - el, ip, co, & + el, ip, co, ma, & h, ce materials => config_material%get('material') @@ -153,6 +161,17 @@ subroutine parse() enddo + allocate(material_orientation0_2(materials%length)) + + do ma = 1, materials%length + material => materials%get(ma) + constituents => material%get('constituents') + allocate(material_orientation0_2(ma)%data(constituents%length)) + do co = 1, constituents%length + call material_orientation0_2(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) + enddo + enddo + end subroutine parse diff --git a/src/phase.f90 b/src/phase.f90 index f0b5c4e22..9d4a423d0 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -19,8 +19,9 @@ module phase implicit none private - type(Rotation), dimension(:,:,:), allocatable :: & - material_orientation0 !< initial orientation of each grain,IP,element + type(tRotationContainer), dimension(:), allocatable :: & + phase_orientation0, & + phase_orientation type(rotation), dimension(:,:,:), allocatable :: & crystallite_orientation !< current orientation @@ -346,7 +347,7 @@ contains subroutine phase_init integer :: & - ph + ph, ce, co, ma class (tNode), pointer :: & debug_constitutive, & materials, & @@ -367,6 +368,25 @@ subroutine phase_init materials => config_material%get('material') phases => config_material%get('phase') + + allocate(phase_orientation0(phases%length)) + do ph = 1,phases%length + allocate(phase_orientation0(ph)%data(count(material_phaseID==ph))) + enddo + + do ce = 1, size(material_phaseID,2) + ma = discretization_materialAt((ce-1)/discretization_nIPs+1) + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + ph = material_phaseID(co,ce) + phase_orientation0(ph)%data(material_phaseEntry(co,ce)) = material_orientation0_2(ma)%data(co) + enddo + enddo + + allocate(phase_orientation(phases%length)) + do ph = 1,phases%length + phase_orientation(ph)%data = phase_orientation0(ph)%data + enddo + call mechanical_init(materials,phases) call damage_init call thermal_init(phases) @@ -602,7 +622,7 @@ function crystallite_push33ToRef(co,ce, tensor33) ph = material_phaseID(co,ce) en = material_phaseEntry(co,ce) - T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct? + T = matmul(phase_orientation0(ph)%data(en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct? crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index cbaf149ec..de23d39fa 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -68,7 +68,7 @@ submodule(phase) mechanical dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient end subroutine phase_hooke_SandItsTangents - + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -204,7 +204,6 @@ module subroutine mechanical_init(materials,phases) num_crystallite, & material, & constituents, & - constituent, & phase, & mech @@ -228,8 +227,6 @@ module subroutine mechanical_init(materials,phases) allocate(phase_mechanical_P(phases%length)) allocate(phase_mechanical_S0(phases%length)) - allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseEntry))) - do ph = 1, phases%length Nmembers = count(material_phaseID == ph) @@ -260,15 +257,11 @@ module subroutine mechanical_init(materials,phases) do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) material => materials%get(discretization_materialAt(el)) - constituents => material%get('constituents') - constituent => constituents%get(co) ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) - call material_orientation0(co,ph,en)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) - - phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = material_orientation0(co,ph,en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) + phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_orientation0(ph)%data(en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) & / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal) phase_mechanical_Fi0(ph)%data(1:3,1:3,en) = math_I3 diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 8eb90f547..03376721f 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1450,8 +1450,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) elseif (prm%chi_GB >= 0.0_pReal) then !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) - if (any(dNeq(material_orientation0(1,ph,en)%asQuaternion(), & - material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. & + if (any(dNeq(phase_orientation0(ph)%data(en)%asQuaternion(), & + phase_orientation0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & (.not. phase_localPlasticity(neighbor_phase))) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) else From d41f3a54900600f09532ac93e35f09804254c986 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 18:42:18 +0200 Subject: [PATCH 08/18] consistent data structure --- src/material.f90 | 11 +- src/phase.f90 | 21 ++-- src/phase_mechanical.f90 | 142 ++++++++++------------ src/phase_mechanical_plastic_nonlocal.f90 | 4 +- 4 files changed, 81 insertions(+), 97 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index c62442fcd..75ddceba5 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -20,7 +20,7 @@ module material type(Rotation), dimension(:), allocatable :: data end type - type(tRotationContainer), dimension(:), allocatable :: material_orientation0_2 + type(tRotationContainer), dimension(:), allocatable :: material_orientation0 integer, dimension(:), allocatable, public, protected :: & homogenization_Nconstituents !< number of grains in each homogenization @@ -46,7 +46,7 @@ module material public :: & tRotationContainer, & - material_orientation0_2, & + material_orientation0, & material_init contains @@ -161,14 +161,15 @@ subroutine parse() enddo - allocate(material_orientation0_2(materials%length)) + allocate(material_orientation0(materials%length)) do ma = 1, materials%length material => materials%get(ma) constituents => material%get('constituents') - allocate(material_orientation0_2(ma)%data(constituents%length)) + allocate(material_orientation0(ma)%data(constituents%length)) do co = 1, constituents%length - call material_orientation0_2(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) + constituent => constituents%get(co) + call material_orientation0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) enddo enddo diff --git a/src/phase.f90 b/src/phase.f90 index 9d4a423d0..9dd8df094 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -23,9 +23,6 @@ module phase phase_orientation0, & phase_orientation - type(rotation), dimension(:,:,:), allocatable :: & - crystallite_orientation !< current orientation - type :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data end type @@ -259,8 +256,7 @@ module phase ph, & i, & e - type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: & - orientation !< crystal orientation + type(tRotationContainer), dimension(:), intent(in) :: orientation end subroutine plastic_nonlocal_updateCompatibility module subroutine plastic_dependentState(co,ip,el) @@ -378,7 +374,7 @@ subroutine phase_init ma = discretization_materialAt((ce-1)/discretization_nIPs+1) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) ph = material_phaseID(co,ce) - phase_orientation0(ph)%data(material_phaseEntry(co,ce)) = material_orientation0_2(ma)%data(co) + phase_orientation0(ph)%data(material_phaseEntry(co,ce)) = material_orientation0(ma)%data(co) enddo enddo @@ -523,8 +519,6 @@ subroutine crystallite_init() iMax = discretization_nIPs eMax = discretization_Nelems - allocate(crystallite_orientation(cMax,iMax,eMax)) - num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal) @@ -594,13 +588,16 @@ subroutine crystallite_orientations(co,ip,el) ip, & !< counter in integration point loop el !< counter in element loop + integer :: ph, en - call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(& - mechanical_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))))) + + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + + call phase_orientation(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en)))) if (plasticState(material_phaseAt(1,el))%nonlocal) & - call plastic_nonlocal_updateCompatibility(crystallite_orientation, & - material_phaseAt(1,el),ip,el) + call plastic_nonlocal_updateCompatibility(phase_orientation,material_phaseAt(1,el),ip,el) end subroutine crystallite_orientations diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index de23d39fa..ed3a737b1 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -248,18 +248,14 @@ module subroutine mechanical_init(materials,phases) phase => phases%get(ph) mech => phase%get('mechanical') #if defined(__GFORTRAN__) - output_constituent(ph)%label = output_as1dString(mech) + output_constituent(ph)%label = output_as1dString(mech) #else - output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray) + output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray) #endif enddo - do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - material => materials%get(discretization_materialAt(el)) - - ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) - en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + do ph = 1, phases%length + do en = 1, count(material_phaseID == ph) phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_orientation0(ph)%data(en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) & @@ -268,13 +264,12 @@ module subroutine mechanical_init(materials,phases) phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3 phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), & - phase_mechanical_Fp0(ph)%data(1:3,1:3,en))) ! assuming that euler angles are given in internal strain free configuration - phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) - phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) - phase_mechanical_F(ph)%data(1:3,1:3,en) = phase_mechanical_F0(ph)%data(1:3,1:3,en) - + phase_mechanical_Fp0(ph)%data(1:3,1:3,en))) ! assuming that euler angles are given in internal strain free configuration enddo - enddo; enddo + phase_mechanical_Fp(ph)%data = phase_mechanical_Fp0(ph)%data + phase_mechanical_Fi(ph)%data = phase_mechanical_Fi0(ph)%data + phase_mechanical_F(ph)%data = phase_mechanical_F0(ph)%data + enddo ! initialize elasticity @@ -903,80 +898,71 @@ subroutine crystallite_results(group,ph) 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_ISO_ID) - structureLabel = 'aP' - case(lattice_FCC_ID) - structureLabel = 'cF' - case(lattice_BCC_ID) - structureLabel = 'cI' - case(lattice_BCT_ID) - structureLabel = 'tI' - case(lattice_HEX_ID) - structureLabel = 'hP' - case(lattice_ORT_ID) - structureLabel = 'oP' - end select - selected_rotations = select_rotations(crystallite_orientation,ph) - 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') + select case(lattice_structure(ph)) + case(lattice_ISO_ID) + structureLabel = 'aP' + case(lattice_FCC_ID) + structureLabel = 'cF' + case(lattice_BCC_ID) + structureLabel = 'cI' + case(lattice_BCT_ID) + structureLabel = 'tI' + case(lattice_HEX_ID) + structureLabel = 'hP' + case(lattice_ORT_ID) + structureLabel = 'oP' + 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 contains !-------------------------------------------------------------------------------------------------- -!> @brief select rotations for output +!> @brief Convert orientation for output: ToDo: implement in HDF5/results !-------------------------------------------------------------------------------------------------- - function select_rotations(dataset,ph) + function select_rotations(dataset) - integer, intent(in) :: ph - type(rotation), dimension(:,:,:), intent(in) :: dataset - real(pReal), dimension(4,count(material_phaseID==ph)) :: select_rotations - integer :: el,ip,co,j + type(rotation), dimension(:), intent(in) :: dataset + real(pReal), dimension(4,size(dataset,1)) :: select_rotations + integer :: en - j=0 - do el = 1, size(material_phaseAt,2) - do ip = 1, discretization_nIPs - do co = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains - if (material_phaseAt(co,el) == ph) then - j = j + 1 - select_rotations(1:4,j) = dataset(co,ip,el)%asQuaternion() - endif - enddo - enddo + do en = 1, size(dataset,1) + select_rotations(:,en) = dataset(en)%asQuaternion() enddo end function select_rotations diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 03376721f..0eccf767e 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1394,7 +1394,7 @@ end function rhoDotFlux !-------------------------------------------------------------------------------------------------- module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) - type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: & + type(tRotationContainer), dimension(:), intent(in) :: & orientation ! crystal orientation integer, intent(in) :: & ph, & @@ -1464,7 +1464,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. - mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) + mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) mySlipSystems: do s1 = 1,ns neighborSlipSystems: do s2 = 1,ns my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & From ee80efd705200d9fba1f0f31f3539981a5ff6ae8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 13:18:33 +0200 Subject: [PATCH 09/18] using new mappings --- src/material.f90 | 4 +- src/phase_damage_anisobrittle.f90 | 2 +- src/results.f90 | 62 ++++++++++++------------------- 3 files changed, 27 insertions(+), 41 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index 75ddceba5..c9f6e8950 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -68,8 +68,8 @@ subroutine material_init(restart) if (.not. restart) then call results_openJobFile - call results_mapping_phase(material_phaseAt,material_phaseMemberAt,material_name_phase) - call results_mapping_homogenization(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization) + call results_mapping_phase(material_phaseID,material_phaseMemberAt,material_name_phase) + call results_mapping_homogenization(material_homogenizationID,material_homogenizationMemberAt,material_name_homogenization) call results_closeJobFile endif diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index ddcca9044..40fe2abe1 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -92,7 +92,7 @@ module function anisobrittle_init() result(mySources) if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' - Nmembers = count(material_phaseAt==p) * discretization_nIPs + Nmembers = count(material_phaseID==p) call phase_allocateState(damageState(p),Nmembers,1,1,0) damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' diff --git a/src/results.f90 b/src/results.f90 index f0af4a1cb..c34867c6f 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -415,16 +415,15 @@ end subroutine results_writeTensorDataset_int !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine results_mapping_phase(phaseAt,memberAtLocal,label) +subroutine results_mapping_phase(phaseID,memberAtLocal,label) - integer, dimension(:,:), intent(in) :: phaseAt !< phase section at (constituent,element) - integer, dimension(:,:,:), intent(in) :: memberAtLocal !< phase member at (constituent,IP,element) + integer, dimension(:,:), intent(in) :: phaseID !< phase ID at (constituent,ce) + integer, dimension(:,:,:), intent(in) :: memberAtLocal !< phase entry at (constituent,IP,element) character(len=*), dimension(:), intent(in) :: label !< label of each phase section integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2),size(memberAtLocal,3)) :: & - phaseAtMaterialpoint, & memberAtGlobal - integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in member counting per process + integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in entry counting per process integer, dimension(0:worldsize-1) :: writeSize !< amount of data written per process integer(HSIZE_T), dimension(2) :: & myShape, & !< shape of the dataset (this process) @@ -451,10 +450,10 @@ subroutine results_mapping_phase(phaseAt,memberAtLocal,label) if(hdferr < 0) error stop 'HDF5 error' memberOffset = 0 do i=1, size(label) - memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAtLocal,2) ! number of points/instance of this process + memberOffset(i,worldrank) = count(phaseID == i) ! number of entries of this process enddo writeSize = 0 - writeSize(worldrank) = size(memberAtLocal(1,:,:)) ! total number of points by this process + writeSize(worldrank) = size(memberAtLocal(1,:,:)) ! total number of entries of this process !-------------------------------------------------------------------------------------------------- ! MPI settings and communication @@ -469,21 +468,15 @@ subroutine results_mapping_phase(phaseAt,memberAtLocal,label) if(ierr /= 0) error stop 'MPI error' #endif - myShape = int([size(phaseAt,1),writeSize(worldrank)], HSIZE_T) - myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T) - totalShape = int([size(phaseAt,1),sum(writeSize)], HSIZE_T) - - -!--------------------------------------------------------------------------------------------------- -! expand phaseAt to consider IPs (is not stored per IP) - do i = 1, size(phaseAtMaterialpoint,2) - phaseAtMaterialpoint(:,i,:) = phaseAt - enddo + myShape = int([size(phaseID,1),writeSize(worldrank)], HSIZE_T) + myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T) + totalShape = int([size(phaseID,1),sum(writeSize)], HSIZE_T) !--------------------------------------------------------------------------------------------------- ! renumber member from my process to all processes do i = 1, size(label) - where(phaseAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1 ! convert to 0-based + where(reshape(phaseID,shape(memberAtGlobal)) == i) & + memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1 ! 0-based enddo !--------------------------------------------------------------------------------------------------- @@ -540,7 +533,7 @@ subroutine results_mapping_phase(phaseAt,memberAtLocal,label) call h5dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' - call h5dwrite_f(dset_id, label_id, reshape(label(pack(phaseAtMaterialpoint,.true.)),myShape), & + call h5dwrite_f(dset_id, label_id, reshape(label(pack(phaseID,.true.)),myShape), & myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if(hdferr < 0) error stop 'HDF5 error' call h5dwrite_f(dset_id, entry_id, reshape(pack(memberAtGlobal,.true.),myShape), & @@ -572,16 +565,15 @@ end subroutine results_mapping_phase !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) +subroutine results_mapping_homogenization(homogenizationID,memberAtLocal,label) - integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element) - integer, dimension(:,:), intent(in) :: memberAtLocal !< homogenization member at (IP,element) + integer, dimension(:), intent(in) :: homogenizationID !< homogenization ID at (cell) + integer, dimension(:,:), intent(in) :: memberAtLocal !< homogenization entry at (IP,element) character(len=*), dimension(:), intent(in) :: label !< label of each homogenization section integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2)) :: & - homogenizationAtMaterialpoint, & memberAtGlobal - integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in member counting per process + integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in entry counting per process integer, dimension(0:worldsize-1) :: writeSize !< amount of data written per process integer(HSIZE_T), dimension(1) :: & myShape, & !< shape of the dataset (this process) @@ -609,10 +601,10 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) if(hdferr < 0) error stop 'HDF5 error' memberOffset = 0 do i=1, size(label) - memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAtLocal,1) ! number of points/instance of this process + memberOffset(i,worldrank) = count(homogenizationID == i) ! number of entries of this process enddo writeSize = 0 - writeSize(worldrank) = size(memberAtLocal) ! total number of points by this process + writeSize(worldrank) = size(memberAtLocal) ! total number of entries of this process !-------------------------------------------------------------------------------------------------- ! MPI settings and communication @@ -627,21 +619,15 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) if(ierr /= 0) error stop 'MPI error' #endif - myShape = int([writeSize(worldrank)], HSIZE_T) - myOffset = int([sum(writeSize(0:worldrank-1))], HSIZE_T) - totalShape = int([sum(writeSize)], HSIZE_T) - - -!--------------------------------------------------------------------------------------------------- -! expand phaseAt to consider IPs (is not stored per IP) - do i = 1, size(homogenizationAtMaterialpoint,1) - homogenizationAtMaterialpoint(i,:) = homogenizationAt - enddo + myShape = int([writeSize(worldrank)], HSIZE_T) + myOffset = int([sum(writeSize(0:worldrank-1))], HSIZE_T) + totalShape = int([sum(writeSize)], HSIZE_T) !--------------------------------------------------------------------------------------------------- ! renumber member from my process to all processes do i = 1, size(label) - where(homogenizationAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1 ! convert to 0-based + where(reshape(homogenizationID,shape(memberAtGlobal)) == i) & + memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1 ! 0-based enddo !--------------------------------------------------------------------------------------------------- @@ -698,7 +684,7 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label) call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' - call h5dwrite_f(dset_id, label_id, reshape(label(pack(homogenizationAtMaterialpoint,.true.)),myShape), & + call h5dwrite_f(dset_id, label_id, reshape(label(pack(homogenizationID,.true.)),myShape), & myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if(hdferr < 0) error stop 'HDF5 error' call h5dwrite_f(dset_id, entry_id, reshape(pack(memberAtGlobal,.true.),myShape), & From af9fa9e9a1b1d658df10ce5aae7d8d8d1ee1dc27 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 22 May 2021 13:36:23 +0200 Subject: [PATCH 10/18] using current naming scheme --- src/material.f90 | 4 +-- src/results.f90 | 63 ++++++++++++++++++++++++------------------------ 2 files changed, 33 insertions(+), 34 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index c9f6e8950..662dc2fa2 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -68,8 +68,8 @@ subroutine material_init(restart) if (.not. restart) then call results_openJobFile - call results_mapping_phase(material_phaseID,material_phaseMemberAt,material_name_phase) - call results_mapping_homogenization(material_homogenizationID,material_homogenizationMemberAt,material_name_homogenization) + call results_mapping_phase(material_phaseID,material_phaseEntry,material_name_phase) + call results_mapping_homogenization(material_homogenizationID,material_homogenizationEntry,material_name_homogenization) call results_closeJobFile endif diff --git a/src/results.f90 b/src/results.f90 index c34867c6f..31b2db856 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -415,15 +415,15 @@ end subroutine results_writeTensorDataset_int !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine results_mapping_phase(phaseID,memberAtLocal,label) +subroutine results_mapping_phase(ID,entry,label) - integer, dimension(:,:), intent(in) :: phaseID !< phase ID at (constituent,ce) - integer, dimension(:,:,:), intent(in) :: memberAtLocal !< phase entry at (constituent,IP,element) - character(len=*), dimension(:), intent(in) :: label !< label of each phase section + integer, dimension(:,:), intent(in) :: ID !< phase ID at (co,ce) + integer, dimension(:,:), intent(in) :: entry !< phase entry at (co,ce) + character(len=*), dimension(:), intent(in) :: label !< label of each phase section - integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2),size(memberAtLocal,3)) :: & - memberAtGlobal - integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in entry counting per process + integer, dimension(size(entry,1),size(entry,2)) :: & + entryGlobal + integer, dimension(size(label),0:worldsize-1) :: entryOffset !< offset in entry counting per process integer, dimension(0:worldsize-1) :: writeSize !< amount of data written per process integer(HSIZE_T), dimension(2) :: & myShape, & !< shape of the dataset (this process) @@ -444,16 +444,17 @@ subroutine results_mapping_phase(phaseID,memberAtLocal,label) integer(SIZE_T) :: type_size_string, type_size_int integer :: hdferr, ierr, i + !-------------------------------------------------------------------------------------------------- ! prepare MPI communication (transparent for non-MPI runs) call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' - memberOffset = 0 + entryOffset = 0 do i=1, size(label) - memberOffset(i,worldrank) = count(phaseID == i) ! number of entries of this process + entryOffset(i,worldrank) = count(ID == i) ! number of entries of this process enddo writeSize = 0 - writeSize(worldrank) = size(memberAtLocal(1,:,:)) ! total number of entries of this process + writeSize(worldrank) = size(entry(1,:)) ! total number of entries of this process !-------------------------------------------------------------------------------------------------- ! MPI settings and communication @@ -464,19 +465,18 @@ subroutine results_mapping_phase(phaseID,memberAtLocal,label) call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process if(ierr /= 0) error stop 'MPI error' - call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process + call MPI_allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process if(ierr /= 0) error stop 'MPI error' #endif - myShape = int([size(phaseID,1),writeSize(worldrank)], HSIZE_T) + myShape = int([size(ID,1),writeSize(worldrank)], HSIZE_T) myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T) - totalShape = int([size(phaseID,1),sum(writeSize)], HSIZE_T) + totalShape = int([size(ID,1),sum(writeSize)], HSIZE_T) !--------------------------------------------------------------------------------------------------- ! renumber member from my process to all processes do i = 1, size(label) - where(reshape(phaseID,shape(memberAtGlobal)) == i) & - memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1 ! 0-based + where(ID == i) entryGlobal = entry + sum(entryOffset(i,0:worldrank-1)) -1 ! 0-based enddo !--------------------------------------------------------------------------------------------------- @@ -533,10 +533,10 @@ subroutine results_mapping_phase(phaseID,memberAtLocal,label) call h5dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' - call h5dwrite_f(dset_id, label_id, reshape(label(pack(phaseID,.true.)),myShape), & + call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), & myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if(hdferr < 0) error stop 'HDF5 error' - call h5dwrite_f(dset_id, entry_id, reshape(pack(memberAtGlobal,.true.),myShape), & + call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), & myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if(hdferr < 0) error stop 'HDF5 error' @@ -565,15 +565,15 @@ end subroutine results_mapping_phase !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine results_mapping_homogenization(homogenizationID,memberAtLocal,label) +subroutine results_mapping_homogenization(ID,entry,label) - integer, dimension(:), intent(in) :: homogenizationID !< homogenization ID at (cell) - integer, dimension(:,:), intent(in) :: memberAtLocal !< homogenization entry at (IP,element) - character(len=*), dimension(:), intent(in) :: label !< label of each homogenization section + integer, dimension(:), intent(in) :: ID !< homogenization ID at (ce) + integer, dimension(:), intent(in) :: entry !< homogenization entry at (ce) + character(len=*), dimension(:), intent(in) :: label !< label of each homogenization section - integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2)) :: & - memberAtGlobal - integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in entry counting per process + integer, dimension(size(entry,1)) :: & + entryGlobal + integer, dimension(size(label),0:worldsize-1) :: entryOffset !< offset in entry counting per process integer, dimension(0:worldsize-1) :: writeSize !< amount of data written per process integer(HSIZE_T), dimension(1) :: & myShape, & !< shape of the dataset (this process) @@ -599,12 +599,12 @@ subroutine results_mapping_homogenization(homogenizationID,memberAtLocal,label) ! prepare MPI communication (transparent for non-MPI runs) call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' - memberOffset = 0 + entryOffset = 0 do i=1, size(label) - memberOffset(i,worldrank) = count(homogenizationID == i) ! number of entries of this process + entryOffset(i,worldrank) = count(ID == i) ! number of entries of this process enddo writeSize = 0 - writeSize(worldrank) = size(memberAtLocal) ! total number of entries of this process + writeSize(worldrank) = size(entry) ! total number of entries of this process !-------------------------------------------------------------------------------------------------- ! MPI settings and communication @@ -615,7 +615,7 @@ subroutine results_mapping_homogenization(homogenizationID,memberAtLocal,label) call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process if(ierr /= 0) error stop 'MPI error' - call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process + call MPI_allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get offset at each process if(ierr /= 0) error stop 'MPI error' #endif @@ -626,8 +626,7 @@ subroutine results_mapping_homogenization(homogenizationID,memberAtLocal,label) !--------------------------------------------------------------------------------------------------- ! renumber member from my process to all processes do i = 1, size(label) - where(reshape(homogenizationID,shape(memberAtGlobal)) == i) & - memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1 ! 0-based + where(ID == i) entryGlobal = entry + sum(entryOffset(i,0:worldrank-1)) - 1 ! 0-based enddo !--------------------------------------------------------------------------------------------------- @@ -684,10 +683,10 @@ subroutine results_mapping_homogenization(homogenizationID,memberAtLocal,label) call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' - call h5dwrite_f(dset_id, label_id, reshape(label(pack(homogenizationID,.true.)),myShape), & + call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), & myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if(hdferr < 0) error stop 'HDF5 error' - call h5dwrite_f(dset_id, entry_id, reshape(pack(memberAtGlobal,.true.),myShape), & + call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), & myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if(hdferr < 0) error stop 'HDF5 error' From 72ab936ec38f8fc63b1c993b067f9ab70bc842d0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 23 May 2021 00:10:46 +0200 Subject: [PATCH 11/18] cleaning --- src/config.f90 | 2 +- src/material.f90 | 30 ++++++++++++------------------ src/math.f90 | 3 +-- src/parallelization.f90 | 6 +++--- src/phase.f90 | 4 ++-- src/prec.f90 | 2 +- src/results.f90 | 3 +-- 7 files changed, 21 insertions(+), 29 deletions(-) diff --git a/src/config.f90 b/src/config.f90 index 02b16f2a2..7d75cb444 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -92,7 +92,7 @@ end subroutine parse_debug !-------------------------------------------------------------------------------------------------- !> @brief Deallocate config_material. -!ToDo: deallocation of numerics debug (optional) +!ToDo: deallocation of numerics and debug (optional) !-------------------------------------------------------------------------------------------------- subroutine config_deallocate diff --git a/src/material.f90 b/src/material.f90 index 662dc2fa2..9d67fc34c 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -32,18 +32,15 @@ module material material_name_homogenization !< name of each homogenization integer, dimension(:), allocatable, public, protected :: & ! (elem) - material_homogenizationAt, & !< homogenization ID of each element - material_homogenizationID, & !< per cell - material_homogenizationEntry !< per cell - integer, dimension(:,:), allocatable :: & ! (ip,elem) - material_homogenizationMemberAt !< position of the element within its homogenization instance + material_homogenizationAt, & !< homogenization ID of each element TODO: remove + material_homogenizationID, & !< per cell TODO: material_ID_homogenization + material_homogenizationEntry !< per cell TODO: material_entry_homogenization integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) - material_phaseAt, & !< phase ID of each element - material_phaseID, & !< per (constituent,cell) - material_phaseEntry !< per (constituent,cell) + material_phaseAt, & !< phase ID of each element TODO: remove + material_phaseID, & !< per (constituent,cell) TODO: material_ID_phase + material_phaseEntry !< per (constituent,cell) TODO: material_entry_phase integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) - material_phaseMemberAt !< position of the element within its phase instance - + material_phaseMemberAt !TODO: remove public :: & tRotationContainer, & material_orientation0, & @@ -118,7 +115,6 @@ subroutine parse() allocate(counterHomogenization(homogenizations%length),source=0) allocate(material_homogenizationAt(discretization_Nelems),source=0) - allocate(material_homogenizationMemberAt(discretization_nIPs,discretization_Nelems),source=0) allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0) allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) @@ -136,9 +132,8 @@ subroutine parse() do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1 - material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) - material_homogenizationID(ce) = material_homogenizationAt(el) - material_homogenizationEntry(ce) = material_homogenizationMemberAt(ip,el) + material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationAt(el)) + material_homogenizationID(ce) = material_homogenizationAt(el) enddo frac = 0.0_pReal @@ -150,10 +145,9 @@ subroutine parse() do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1 - material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) - - material_phaseID(co,ce) = material_phaseAt(co,el) - material_phaseEntry(co,ce) = material_phaseMemberAt(co,ip,el) + material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) + material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el)) + material_phaseID(co,ce) = material_phaseAt(co,el) enddo enddo diff --git a/src/math.f90 b/src/math.f90 index 6ef942677..5b2cedacb 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1140,8 +1140,7 @@ end function math_areaTriangle !-------------------------------------------------------------------------------------------------- -!> @brief limits a scalar value to a certain range (either one or two sided) -! Will return NaN if left > right +!> @brief Limit a scalar value to a certain range (either one or two sided). !-------------------------------------------------------------------------------------------------- real(pReal) pure elemental function math_clip(a, left, right) diff --git a/src/parallelization.f90 b/src/parallelization.f90 index e850ff04d..ee8e71006 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -27,7 +27,7 @@ module parallelization contains !-------------------------------------------------------------------------------------------------- -!> @brief calls subroutines that reads material, numerics and debug configuration files +!> @brief Initialize shared memory (openMP) and distributed memory (MPI) parallelization. !-------------------------------------------------------------------------------------------------- subroutine parallelization_init @@ -42,8 +42,8 @@ subroutine parallelization_init ! If openMP is enabled, check if the MPI libary supports it and initialize accordingly. ! Otherwise, the first call to PETSc will do the initialization. call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,err) - if (err /= 0) error stop 'MPI init failed' - if (threadLevel @brief report precision and do self test +!> @brief Report precision and do self test. !-------------------------------------------------------------------------------------------------- subroutine prec_init diff --git a/src/results.f90 b/src/results.f90 index 31b2db856..97a3b2489 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -423,7 +423,7 @@ subroutine results_mapping_phase(ID,entry,label) integer, dimension(size(entry,1),size(entry,2)) :: & entryGlobal - integer, dimension(size(label),0:worldsize-1) :: entryOffset !< offset in entry counting per process + integer, dimension(size(label),0:worldsize-1) :: entryOffset !< offset in entry counting per process integer, dimension(0:worldsize-1) :: writeSize !< amount of data written per process integer(HSIZE_T), dimension(2) :: & myShape, & !< shape of the dataset (this process) @@ -718,7 +718,6 @@ end subroutine results_mapping_homogenization !-------------------------------------------------------------------------------------------------- subroutine executionStamp(path,description,SIunit) - character(len=*), intent(in) :: path,description character(len=*), intent(in), optional :: SIunit From aa9e4aa841f93f1333969457e3dae47d6644e21b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 23 May 2021 00:19:46 +0200 Subject: [PATCH 12/18] not needed --- examples/mesh/numerics.config | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 examples/mesh/numerics.config diff --git a/examples/mesh/numerics.config b/examples/mesh/numerics.config deleted file mode 100644 index 330114825..000000000 --- a/examples/mesh/numerics.config +++ /dev/null @@ -1,3 +0,0 @@ -residualStiffness 0.001 -charLength 0.02 -petsc_options -mech_snes_type newtonls -mech_ksp_type fgmres -mech_pc_type ml -mech_ksp_monitor \ No newline at end of file From 8dcf4354e1d9fcd6bee8fd5d917c017877331723 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 23 May 2021 08:19:43 +0200 Subject: [PATCH 13/18] more logical structure --- src/results.f90 | 64 ++++++++++++++++++++++++------------------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/src/results.f90 b/src/results.f90 index 97a3b2489..29f1a5fd7 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -445,20 +445,22 @@ subroutine results_mapping_phase(ID,entry,label) integer :: hdferr, ierr, i -!-------------------------------------------------------------------------------------------------- -! prepare MPI communication (transparent for non-MPI runs) - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - entryOffset = 0 - do i=1, size(label) - entryOffset(i,worldrank) = count(ID == i) ! number of entries of this process - enddo + entryGlobal = entry -1 ! 0-based + writeSize = 0 writeSize(worldrank) = size(entry(1,:)) ! total number of entries of this process + entryOffset = 0 + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + !-------------------------------------------------------------------------------------------------- ! MPI settings and communication #ifdef PETSc + do i=1, size(label) + entryOffset(i,worldrank) = count(ID == i) ! number of entries of this process + enddo + call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) if(hdferr < 0) error stop 'HDF5 error' @@ -467,18 +469,16 @@ subroutine results_mapping_phase(ID,entry,label) call MPI_allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process if(ierr /= 0) error stop 'MPI error' + + do i = 1, size(label) + where(ID == i) entryGlobal = entryGlobal + sum(entryOffset(i,0:worldrank-1)) + enddo #endif myShape = int([size(ID,1),writeSize(worldrank)], HSIZE_T) myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T) totalShape = int([size(ID,1),sum(writeSize)], HSIZE_T) -!--------------------------------------------------------------------------------------------------- -! renumber member from my process to all processes - do i = 1, size(label) - where(ID == i) entryGlobal = entry + sum(entryOffset(i,0:worldrank-1)) -1 ! 0-based - enddo - !--------------------------------------------------------------------------------------------------- ! compound type: name of phase section + position/index within results array call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) @@ -595,20 +595,22 @@ subroutine results_mapping_homogenization(ID,entry,label) integer :: hdferr, ierr, i -!-------------------------------------------------------------------------------------------------- -! prepare MPI communication (transparent for non-MPI runs) - call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if(hdferr < 0) error stop 'HDF5 error' - entryOffset = 0 - do i=1, size(label) - entryOffset(i,worldrank) = count(ID == i) ! number of entries of this process - enddo + entryGlobal = entry -1 ! 0-based + writeSize = 0 writeSize(worldrank) = size(entry) ! total number of entries of this process + entryOffset = 0 + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + !-------------------------------------------------------------------------------------------------- ! MPI settings and communication #ifdef PETSc + do i=1, size(label) + entryOffset(i,worldrank) = count(ID == i) ! number of entries of this process + enddo + call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) if(hdferr < 0) error stop 'HDF5 error' @@ -617,18 +619,16 @@ subroutine results_mapping_homogenization(ID,entry,label) call MPI_allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get offset at each process if(ierr /= 0) error stop 'MPI error' + + do i = 1, size(label) + where(ID == i) entryGlobal = entryGlobal + sum(entryOffset(i,0:worldrank-1)) + enddo #endif myShape = int([writeSize(worldrank)], HSIZE_T) myOffset = int([sum(writeSize(0:worldrank-1))], HSIZE_T) totalShape = int([sum(writeSize)], HSIZE_T) -!--------------------------------------------------------------------------------------------------- -! renumber member from my process to all processes - do i = 1, size(label) - where(ID == i) entryGlobal = entry + sum(entryOffset(i,0:worldrank-1)) - 1 ! 0-based - enddo - !--------------------------------------------------------------------------------------------------- ! compound type: name of phase section + position/index within results array call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) @@ -722,14 +722,14 @@ subroutine executionStamp(path,description,SIunit) character(len=*), intent(in), optional :: SIunit - if (HDF5_objectExists(resultsFile,path)) & - call HDF5_addAttribute(resultsFile,'description',description,path) - if (HDF5_objectExists(resultsFile,path) .and. present(SIunit)) & - call HDF5_addAttribute(resultsFile,'unit',SIunit,path) if (HDF5_objectExists(resultsFile,path)) & call HDF5_addAttribute(resultsFile,'creator','DAMASK '//DAMASKVERSION,path) if (HDF5_objectExists(resultsFile,path)) & call HDF5_addAttribute(resultsFile,'created',now(),path) + if (HDF5_objectExists(resultsFile,path)) & + call HDF5_addAttribute(resultsFile,'description',description,path) + if (HDF5_objectExists(resultsFile,path) .and. present(SIunit)) & + call HDF5_addAttribute(resultsFile,'unit',SIunit,path) end subroutine executionStamp From c4d19691502507cce8864c315896a486251cc20a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 23 May 2021 09:02:57 +0200 Subject: [PATCH 14/18] easier to read and faster for multiphase materials --- src/results.f90 | 38 ++++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/src/results.f90 b/src/results.f90 index 29f1a5fd7..3ab23b163 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -442,7 +442,7 @@ subroutine results_mapping_phase(ID,entry,label) dt_id integer(SIZE_T) :: type_size_string, type_size_int - integer :: hdferr, ierr, i + integer :: hdferr, ierr, ce, co entryGlobal = entry -1 ! 0-based @@ -450,28 +450,31 @@ subroutine results_mapping_phase(ID,entry,label) writeSize = 0 writeSize(worldrank) = size(entry(1,:)) ! total number of entries of this process - entryOffset = 0 call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' !-------------------------------------------------------------------------------------------------- ! MPI settings and communication #ifdef PETSc - do i=1, size(label) - entryOffset(i,worldrank) = count(ID == i) ! number of entries of this process - enddo - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) if(hdferr < 0) error stop 'HDF5 error' call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process if(ierr /= 0) error stop 'MPI error' + entryOffset = 0 + do co = 1, size(ID,1) + do ce = 1, size(ID,2) + entryOffset(ID(co,ce),worldrank) = entryOffset(ID(co,ce),worldrank) +1 + enddo + enddo call MPI_allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process if(ierr /= 0) error stop 'MPI error' - - do i = 1, size(label) - where(ID == i) entryGlobal = entryGlobal + sum(entryOffset(i,0:worldrank-1)) + entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2) + do co = 1, size(ID,1) + do ce = 1, size(ID,2) + entryGlobal(co,ce) = entryGlobal(co,ce) + entryOffset(ID(co,ce),worldrank) + enddo enddo #endif @@ -592,7 +595,7 @@ subroutine results_mapping_homogenization(ID,entry,label) dt_id integer(SIZE_T) :: type_size_string, type_size_int - integer :: hdferr, ierr, i + integer :: hdferr, ierr, i, ce entryGlobal = entry -1 ! 0-based @@ -600,28 +603,27 @@ subroutine results_mapping_homogenization(ID,entry,label) writeSize = 0 writeSize(worldrank) = size(entry) ! total number of entries of this process - entryOffset = 0 call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' !-------------------------------------------------------------------------------------------------- ! MPI settings and communication #ifdef PETSc - do i=1, size(label) - entryOffset(i,worldrank) = count(ID == i) ! number of entries of this process - enddo - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) if(hdferr < 0) error stop 'HDF5 error' call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process if(ierr /= 0) error stop 'MPI error' + entryOffset = 0 + do ce = 1, size(ID,1) + entryOffset(ID(ce),worldrank) = entryOffset(ID(ce),worldrank) +1 + enddo call MPI_allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get offset at each process if(ierr /= 0) error stop 'MPI error' - - do i = 1, size(label) - where(ID == i) entryGlobal = entryGlobal + sum(entryOffset(i,0:worldrank-1)) + entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2) + do ce = 1, size(ID,1) + entryGlobal(ce) = entryGlobal(ce) + entryOffset(ID(ce),worldrank) enddo #endif From 4713e0e85db0bef93c61e21567a6531f67a4c9cc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 23 May 2021 10:10:25 +0200 Subject: [PATCH 15/18] new names and mappings --- src/CPFEM.f90 | 5 ----- src/homogenization.f90 | 11 ++++++----- src/homogenization_damage.f90 | 10 +++++----- src/homogenization_mechanical_RGC.f90 | 20 ++++++++++---------- src/homogenization_mechanical_isostrain.f90 | 8 ++++---- src/homogenization_mechanical_pass.f90 | 8 ++++---- src/material.f90 | 11 ++++------- src/phase.f90 | 4 +++- src/phase_mechanical.f90 | 4 ++-- src/phase_mechanical_plastic_nonlocal.f90 | 1 - 10 files changed, 38 insertions(+), 44 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index d8f4a0c3c..213b56a54 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -177,11 +177,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward - !chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) - ! case (THERMAL_conduction_ID) chosenThermal1 - ! temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & - ! temperature_inp - !end select chosenThermal1 homogenization_F0(1:3,1:3,ce) = ffn homogenization_F(1:3,1:3,ce) = ffn1 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 0d62f83c7..f25a749c2 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -245,11 +245,12 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE !$OMP PARALLEL !$OMP DO PRIVATE(ce,en,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) - ho = material_homogenizationAt(el) - myNgrains = homogenization_Nconstituents(ho) + do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip en = material_homogenizationEntry(ce) + ho = material_homogenizationID(ce) + myNgrains = homogenization_Nconstituents(ho) call phase_restore(ce,.false.) ! wrong name (is more a forward function) @@ -290,12 +291,12 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE !$OMP DO PRIVATE(ho,ph,ce) do el = FEsolving_execElem(1),FEsolving_execElem(2) if (terminallyIll) continue - ho = material_homogenizationAt(el) do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip + ho = material_homogenizationID(ce) call thermal_partition(ce) do co = 1, homogenization_Nconstituents(ho) - ph = material_phaseAt(co,el) + ph = material_phaseID(co,ce) if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then if (.not. terminallyIll) & ! so first signals terminally ill... print*, ' Integration point ', ip,' at element ', el, ' terminally ill' @@ -308,9 +309,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE !$OMP DO PRIVATE(ho,ce) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) - ho = material_homogenizationAt(el) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip + ho = material_homogenizationID(ce) do co = 1, homogenization_Nconstituents(ho) call crystallite_orientations(co,ip,el) enddo diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 0546834fd..b33616a7c 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -39,7 +39,7 @@ module subroutine damage_init() configHomogenization, & configHomogenizationDamage, & num_generic - integer :: ho,Nmaterialpoints + integer :: ho,Nmembers print'(/,a)', ' <<<+- homogenization:damage init -+>>>' @@ -50,7 +50,8 @@ module subroutine damage_init() allocate(current(configHomogenizations%length)) do ho = 1, configHomogenizations%length - allocate(current(ho)%phi(count(material_homogenizationID==ho)), source=1.0_pReal) + Nmembers = count(material_homogenizationID == ho) + allocate(current(ho)%phi(Nmembers), source=1.0_pReal) configHomogenization => configHomogenizations%get(ho) associate(prm => param(ho)) if (configHomogenization%contains('damage')) then @@ -60,10 +61,9 @@ module subroutine damage_init() #else prm%output = configHomogenizationDamage%get_as1dString('output',defaultVal=emptyStringArray) #endif - Nmaterialpoints = count(material_homogenizationAt == ho) damageState_h(ho)%sizeState = 1 - allocate(damageState_h(ho)%state0(1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState_h(ho)%state (1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState_h(ho)%state0(1,Nmembers), source=1.0_pReal) + allocate(damageState_h(ho)%state (1,Nmembers), source=1.0_pReal) else prm%output = emptyStringArray endif diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 745c266d4..6392e00db 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -78,7 +78,7 @@ module subroutine RGC_init(num_homogMech) integer :: & ho, & - Nmaterialpoints, & + Nmembers, & sizeState, nIntFaceTot class (tNode), pointer :: & @@ -161,28 +161,28 @@ module subroutine RGC_init(num_homogMech) prm%D_alpha = homogMech%get_as1dFloat('D_alpha', requiredSize=3) prm%a_g = homogMech%get_as1dFloat('a_g', requiredSize=3) - Nmaterialpoints = count(material_homogenizationAt == ho) + Nmembers = count(material_homogenizationID == ho) nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) & + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) & + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1)) sizeState = nIntFaceTot homogState(ho)%sizeState = sizeState - allocate(homogState(ho)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal) - allocate(homogState(ho)%state (sizeState,Nmaterialpoints), source=0.0_pReal) + allocate(homogState(ho)%state0 (sizeState,Nmembers), source=0.0_pReal) + allocate(homogState(ho)%state (sizeState,Nmembers), source=0.0_pReal) stt%relaxationVector => homogState(ho)%state(1:nIntFaceTot,:) st0%relaxationVector => homogState(ho)%state0(1:nIntFaceTot,:) - allocate(dst%volumeDiscrepancy( Nmaterialpoints), source=0.0_pReal) - allocate(dst%relaxationRate_avg( Nmaterialpoints), source=0.0_pReal) - allocate(dst%relaxationRate_max( Nmaterialpoints), source=0.0_pReal) - allocate(dst%mismatch( 3,Nmaterialpoints), source=0.0_pReal) + allocate(dst%volumeDiscrepancy( Nmembers), source=0.0_pReal) + allocate(dst%relaxationRate_avg( Nmembers), source=0.0_pReal) + allocate(dst%relaxationRate_max( Nmembers), source=0.0_pReal) + allocate(dst%mismatch( 3,Nmembers), source=0.0_pReal) !-------------------------------------------------------------------------------------------------- ! assigning cluster orientations - dependentState(ho)%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints) - !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints) ifort version 18.0.1 crashes (for whatever reason) + dependentState(ho)%orientation = spread(eu2om(prm%a_g*inRad),3,Nmembers) + !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmembers) ifort version 18.0.1 crashes (for whatever reason) end associate diff --git a/src/homogenization_mechanical_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90 index 7b114d04f..1dfac425c 100644 --- a/src/homogenization_mechanical_isostrain.f90 +++ b/src/homogenization_mechanical_isostrain.f90 @@ -15,7 +15,7 @@ module subroutine isostrain_init integer :: & ho, & - Nmaterialpoints + Nmembers print'(/,a)', ' <<<+- homogenization:mechanical:isostrain init -+>>>' @@ -25,10 +25,10 @@ module subroutine isostrain_init do ho = 1, size(homogenization_type) if (homogenization_type(ho) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle - Nmaterialpoints = count(material_homogenizationAt == ho) + Nmembers = count(material_homogenizationID == ho) homogState(ho)%sizeState = 0 - allocate(homogState(ho)%state0(0,Nmaterialpoints)) - allocate(homogState(ho)%state (0,Nmaterialpoints)) + allocate(homogState(ho)%state0(0,Nmembers)) + allocate(homogState(ho)%state (0,Nmembers)) enddo diff --git a/src/homogenization_mechanical_pass.f90 b/src/homogenization_mechanical_pass.f90 index e2e44658a..0728d8a06 100644 --- a/src/homogenization_mechanical_pass.f90 +++ b/src/homogenization_mechanical_pass.f90 @@ -15,7 +15,7 @@ module subroutine pass_init integer :: & ho, & - Nmaterialpoints + Nmembers print'(/,a)', ' <<<+- homogenization:mechanical:pass init -+>>>' @@ -28,10 +28,10 @@ module subroutine pass_init if(homogenization_Nconstituents(ho) /= 1) & call IO_error(211,ext_msg='N_constituents (pass)') - Nmaterialpoints = count(material_homogenizationAt == ho) + Nmembers = count(material_homogenizationID == ho) homogState(ho)%sizeState = 0 - allocate(homogState(ho)%state0(0,Nmaterialpoints)) - allocate(homogState(ho)%state (0,Nmaterialpoints)) + allocate(homogState(ho)%state0(0,Nmembers)) + allocate(homogState(ho)%state (0,Nmembers)) enddo diff --git a/src/material.f90 b/src/material.f90 index 9d67fc34c..7aee676e6 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -32,7 +32,6 @@ module material material_name_homogenization !< name of each homogenization integer, dimension(:), allocatable, public, protected :: & ! (elem) - material_homogenizationAt, & !< homogenization ID of each element TODO: remove material_homogenizationID, & !< per cell TODO: material_ID_homogenization material_homogenizationEntry !< per cell TODO: material_entry_homogenization integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) @@ -40,7 +39,7 @@ module material material_phaseID, & !< per (constituent,cell) TODO: material_ID_phase material_phaseEntry !< per (constituent,cell) TODO: material_entry_phase integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) - material_phaseMemberAt !TODO: remove + material_phaseMemberAt !TODO: remove public :: & tRotationContainer, & material_orientation0, & @@ -114,7 +113,6 @@ subroutine parse() allocate(counterPhase(phases%length),source=0) allocate(counterHomogenization(homogenizations%length),source=0) - allocate(material_homogenizationAt(discretization_Nelems),source=0) allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0) allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) @@ -128,12 +126,11 @@ subroutine parse() material => materials%get(discretization_materialAt(el)) constituents => material%get('constituents') - material_homogenizationAt(el) = homogenizations%getIndex(material%get_asString('homogenization')) do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip - counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1 - material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationAt(el)) - material_homogenizationID(ce) = material_homogenizationAt(el) + material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization')) + counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1 + material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce)) enddo frac = 0.0_pReal diff --git a/src/phase.f90 b/src/phase.f90 index e2cb87332..a72d72d1a 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -501,6 +501,7 @@ subroutine crystallite_init() integer :: & ph, & + ce, & co, & !< counter in integration point component loop ip, & !< counter in integration point loop el, & !< counter in element loop @@ -566,7 +567,8 @@ subroutine crystallite_init() !$OMP PARALLEL DO do el = 1, eMax do ip = 1, iMax - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) + ce = (el-1)*discretization_nIPs + ip + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) call crystallite_orientations(co,ip,el) call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states enddo diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index ed3a737b1..8d07b73df 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1020,8 +1020,8 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), dimension(:), allocatable :: subState0 - ph = material_phaseAt(co,el) - en = material_phaseMemberAt(co,ip,el) + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) sizeDotState = plasticState(ph)%sizeDotState subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 0eccf767e..bde14756e 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -506,7 +506,6 @@ module function plastic_nonlocal_init() result(myPlasticity) end associate if (Nmembers > 0) call stateInit(ini,ph,Nmembers) - plasticState(ph)%state0 = plasticState(ph)%state !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range From 7afb1d28ccdaf31d60a539df30813385cd057840 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 23 May 2021 15:04:55 +0200 Subject: [PATCH 16/18] easier to read --- src/results.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/results.f90 b/src/results.f90 index 3ab23b163..5ba48e995 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -445,17 +445,17 @@ subroutine results_mapping_phase(ID,entry,label) integer :: hdferr, ierr, ce, co - entryGlobal = entry -1 ! 0-based - writeSize = 0 writeSize(worldrank) = size(entry(1,:)) ! total number of entries of this process call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' +#ifndef PETSc + entryGlobal = entry -1 ! 0-based +#else !-------------------------------------------------------------------------------------------------- ! MPI settings and communication -#ifdef PETSc call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) if(hdferr < 0) error stop 'HDF5 error' @@ -473,7 +473,7 @@ subroutine results_mapping_phase(ID,entry,label) entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2) do co = 1, size(ID,1) do ce = 1, size(ID,2) - entryGlobal(co,ce) = entryGlobal(co,ce) + entryOffset(ID(co,ce),worldrank) + entryGlobal(co,ce) = entry(co,ce) -1 + entryOffset(ID(co,ce),worldrank) enddo enddo #endif @@ -483,7 +483,7 @@ subroutine results_mapping_phase(ID,entry,label) totalShape = int([size(ID,1),sum(writeSize)], HSIZE_T) !--------------------------------------------------------------------------------------------------- -! compound type: name of phase section + position/index within results array +! compound type: label(ID) + entry call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr) @@ -598,17 +598,17 @@ subroutine results_mapping_homogenization(ID,entry,label) integer :: hdferr, ierr, i, ce - entryGlobal = entry -1 ! 0-based - writeSize = 0 writeSize(worldrank) = size(entry) ! total number of entries of this process call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' +#ifndef PETSc + entryGlobal = entry -1 ! 0-based +#else !-------------------------------------------------------------------------------------------------- ! MPI settings and communication -#ifdef PETSc call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) if(hdferr < 0) error stop 'HDF5 error' @@ -623,7 +623,7 @@ subroutine results_mapping_homogenization(ID,entry,label) if(ierr /= 0) error stop 'MPI error' entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2) do ce = 1, size(ID,1) - entryGlobal(ce) = entryGlobal(ce) + entryOffset(ID(ce),worldrank) + entryGlobal(ce) = entry(ce) -1 + entryOffset(ID(ce),worldrank) enddo #endif @@ -632,7 +632,7 @@ subroutine results_mapping_homogenization(ID,entry,label) totalShape = int([sum(writeSize)], HSIZE_T) !--------------------------------------------------------------------------------------------------- -! compound type: name of phase section + position/index within results array +! compound type: label(ID) + entry call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr) From 78799f9794cf1e7cf74ad5ad1b87defb510ecf05 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 23 May 2021 18:44:24 +0200 Subject: [PATCH 17/18] fix for openMP --- src/homogenization.f90 | 7 +++---- src/phase.f90 | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index f25a749c2..3162daada 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -236,21 +236,20 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE NiterationMPstate, & ip, & !< integration point number el, & !< element number - myNgrains, co, ce, ho, en, ph + co, ce, ho, en, ph logical :: & converged logical, dimension(2) :: & doneAndHappy !$OMP PARALLEL - !$OMP DO PRIVATE(ce,en,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) + !$OMP DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip en = material_homogenizationEntry(ce) ho = material_homogenizationID(ce) - myNgrains = homogenization_Nconstituents(ho) call phase_restore(ce,.false.) ! wrong name (is more a forward function) @@ -267,7 +266,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) converged = .true. - do co = 1, myNgrains + do co = 1, homogenization_Nconstituents(ho) converged = converged .and. crystallite_stress(dt,co,ip,el) enddo diff --git a/src/phase.f90 b/src/phase.f90 index a72d72d1a..58776229f 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -564,7 +564,7 @@ subroutine crystallite_init() flush(IO_STDOUT) - !$OMP PARALLEL DO + !$OMP PARALLEL DO PRIVATE(ce) do el = 1, eMax do ip = 1, iMax ce = (el-1)*discretization_nIPs + ip From 77a0d5999bd561d3aab346f2f5f4f1889b900c9a Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 28 May 2021 15:03:50 +0200 Subject: [PATCH 18/18] not needed --- src/phase_mechanical.f90 | 5 ----- src/results.f90 | 2 +- 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 88c5d7f0b..57d1098d7 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -193,13 +193,8 @@ module subroutine mechanical_init(materials,phases) phases integer :: & - el, & - ip, & - co, & - ce, & ph, & en, & - stiffDegradationCtr, & Nmembers class(tNode), pointer :: & num_crystallite, & diff --git a/src/results.f90 b/src/results.f90 index 5ba48e995..3c6175b27 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -595,7 +595,7 @@ subroutine results_mapping_homogenization(ID,entry,label) dt_id integer(SIZE_T) :: type_size_string, type_size_int - integer :: hdferr, ierr, i, ce + integer :: hdferr, ierr, ce writeSize = 0