diff --git a/CMakeLists.txt b/CMakeLists.txt index a8e551373..787117524 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,7 +29,7 @@ else() endif() add_definitions("-D${DAMASK_SOLVER}") -set(CMAKE_Fortran_PREPROCESS "ON") +set(CMAKE_Fortran_PREPROCESS "ON") # works only for CMake >= 3.18 # EXPERIMENTAL: This might help to detect HDF5 and FFTW3 in the future if PETSc is not aware of them set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/externalpackages:$ENV{PKG_CONFIG_PATH}") diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index 61ec2143b..43ed64af2 100644 --- a/cmake/Compiler-GNU.cmake +++ b/cmake/Compiler-GNU.cmake @@ -19,6 +19,8 @@ set (STANDARD_CHECK "-std=f2018 -pedantic-errors" ) #------------------------------------------------------------------------------------------------ # Fine tuning compilation options +set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") # preprocessor, needed for CMake < 3.18 + set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIE") # position independent code diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 59ba6d74d..d5e2dff7d 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -22,6 +22,8 @@ set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel") #------------------------------------------------------------------------------------------------ # Fine tuning compilation options +set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") # preprocessor, needed for CMake < 3.18 + set (COMPILE_FLAGS "${COMPILE_FLAGS} -no-ftz") # disable flush underflow to zero, will be set if -O[1,2,3] diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake index 3749b925f..079574c4b 100644 --- a/cmake/Compiler-IntelLLVM.cmake +++ b/cmake/Compiler-IntelLLVM.cmake @@ -24,6 +24,8 @@ set (LINKER_FLAGS "${LINKER_FLAGS} -shared-intel -fc=ifx") #------------------------------------------------------------------------------------------------ # Fine tuning compilation options +set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") # preprocessor, needed for CMake < 3.18 + set (COMPILE_FLAGS "${COMPILE_FLAGS} -no-ftz") # disable flush underflow to zero, will be set if -O[1,2,3] diff --git a/cmake/Compiler-LLVMFlang.cmake b/cmake/Compiler-LLVMFlang.cmake index b28df7027..f6be61b45 100644 --- a/cmake/Compiler-LLVMFlang.cmake +++ b/cmake/Compiler-LLVMFlang.cmake @@ -9,4 +9,4 @@ set (STANDARD_CHECK "-std=f2018 -pedantic" ) #------------------------------------------------------------------------------------------------ # Fine tuning compilation options - +set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") # preprocessor, needed for CMake < 3.18 diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index c8db7d8df..3be9032f3 100644 --- a/src/Marc/materialpoint_Marc.f90 +++ b/src/Marc/materialpoint_Marc.f90 @@ -156,9 +156,11 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, materialpoint_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) else validCalculation - call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip) + call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip, & + (elCP-1)*discretization_nIPs + ip) if (.not. terminallyIll) & - call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP]) + call homogenization_mechanical_response2(dt,(elCP-1)*discretization_nIPs + ip, & + (elCP-1)*discretization_nIPs + ip) terminalIllness: if (terminallyIll) then diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 index 0ea269319..dae4cd4ee 100644 --- a/src/grid/grid_mech_utilities.f90 +++ b/src/grid/grid_mech_utilities.f90 @@ -139,7 +139,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& if (.not. terminallyIll) & call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3) if (.not. terminallyIll) & - call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) + call homogenization_mechanical_response2(Delta_t,1,product(cells(1:2))*cells3) P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 366e38053..726c2fc2e 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -273,6 +273,7 @@ subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end) real(pREAL), intent(in) :: Delta_t !< time increment integer, intent(in) :: & cell_start, cell_end + integer :: & co, ce, ho @@ -296,37 +297,33 @@ end subroutine homogenization_thermal_response !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_mechanical_response2(Delta_t,cell_start,cell_end) real(pREAL), intent(in) :: Delta_t !< time increment - integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + integer, intent(in) :: & + cell_start, cell_end + integer :: & - ip, & !< integration point number - el, & !< element number co, ce, ho - !$OMP PARALLEL DO PRIVATE(ho,ce) - elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) - IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip + !$OMP PARALLEL DO PRIVATE(ho) + do ce = cell_start, cell_end ho = material_ID_homogenization(ce) do co = 1, homogenization_Nconstituents(ho) - call crystallite_orientations(co,ip,el) + call crystallite_orientations(co,ce) end do call mechanical_homogenize(Delta_t,ce) - end do IpLooping3 - end do elementLooping3 + end do !$OMP END PARALLEL DO - end subroutine homogenization_mechanical_response2 !-------------------------------------------------------------------------------------------------- !> @brief writes homogenization results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_result +subroutine homogenization_result() integer :: ho character(len=:), allocatable :: group_base,group @@ -361,7 +358,7 @@ end subroutine homogenization_result !> @brief Forward data after successful increment. ! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- -subroutine homogenization_forward +subroutine homogenization_forward() integer :: ho diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 15a2168a8..0764f3443 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -148,7 +148,7 @@ subroutine utilities_constitutiveResponse(Delta_t,P_av,forwardData) call homogenization_mechanical_response(Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field if (.not. terminallyIll) & - call homogenization_mechanical_response2(Delta_t,[1,mesh_maxNips],[1,mesh_NcpElems]) + call homogenization_mechanical_response2(Delta_t,1,mesh_maxNips*mesh_NcpElems) cutBack = .false. P_av = sum(homogenization_P,dim=3) * wgt diff --git a/src/phase.f90 b/src/phase.f90 index ffa606cc2..7cebda5b7 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -326,11 +326,8 @@ module phase real(pREAL) :: f end function phase_f_T - module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) - integer, intent(in) :: & - ph, & - ip, & - el + module subroutine plastic_nonlocal_updateCompatibility(orientation,ce) + integer, intent(in) :: ce type(tRotationContainer), dimension(:), intent(in) :: orientation end subroutine plastic_nonlocal_updateCompatibility @@ -387,7 +384,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief Initialize constitutive models for individual physics !-------------------------------------------------------------------------------------------------- -subroutine phase_init +subroutine phase_init() integer :: & ph, ce, co, ma @@ -544,25 +541,16 @@ subroutine crystallite_init() integer :: & ce, & co, & !< counter in integration point component loop - ip, & !< counter in integration point loop - el, & !< counter in element loop en, ph - type(tDict), pointer :: & - num_phase, & - phases - phases => config_material%get_dict('phase') - !$OMP PARALLEL DO PRIVATE(ce,ph,en) - do el = 1, discretization_Nelems - do ip = 1, discretization_nIPs - ce = (el-1)*discretization_nIPs + ip - do co = 1,homogenization_Nconstituents(material_ID_homogenization(ce)) - en = material_entry_phase(co,ce) - ph = material_ID_phase(co,ce) - call crystallite_orientations(co,ip,el) - call plastic_dependentState(ph,en) ! update dependent state variables to be consistent with basic states - end do + !$OMP PARALLEL DO PRIVATE(ph,en) + do ce = 1, size(material_ID_homogenization) + do co = 1,homogenization_Nconstituents(material_ID_homogenization(ce)) + ph = material_ID_phase(co,ce) + en = material_entry_phase(co,ce) + call crystallite_orientations(co,ce) + call plastic_dependentState(ph,en) ! update dependent state variables to be consistent with basic states end do end do !$OMP END PARALLEL DO @@ -572,32 +560,30 @@ end subroutine crystallite_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates orientations +!> @brief Update orientations and, if needed, compatibility. !-------------------------------------------------------------------------------------------------- -subroutine crystallite_orientations(co,ip,el) +subroutine crystallite_orientations(co,ce) integer, intent(in) :: & - co, & !< counter in integration point component loop - ip, & !< counter in integration point loop - el !< counter in element loop + co, & + ce integer :: ph, en - ph = material_ID_phase(co,(el-1)*discretization_nIPs + ip) - en = material_entry_phase(co,(el-1)*discretization_nIPs + ip) + ph = material_ID_phase(co,ce) + en = material_entry_phase(co,ce) call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en)))) - if (plasticState(material_ID_phase(1,(el-1)*discretization_nIPs + ip))%nonlocal) & - call plastic_nonlocal_updateCompatibility(phase_O,material_ID_phase(1,(el-1)*discretization_nIPs + ip),ip,el) + if (plasticState(material_ID_phase(1,ce))%nonlocal) call plastic_nonlocal_updateCompatibility(phase_O,ce) end subroutine crystallite_orientations !-------------------------------------------------------------------------------------------------- -!> @brief Map 2nd order tensor to reference config +!> @brief Map 2nd order tensor to reference configuration. !-------------------------------------------------------------------------------------------------- function crystallite_push33ToRef(co,ce, tensor33) @@ -621,15 +607,17 @@ end function crystallite_push33ToRef !-------------------------------------------------------------------------------------------------- -!> @brief determines whether a point is converged +!> @brief Determine whether a point is converged. !-------------------------------------------------------------------------------------------------- logical pure function converged(residuum,state,atol) - real(pREAL), intent(in), dimension(:) ::& + real(pREAL), intent(in), dimension(:) :: & residuum, state, atol + real(pREAL) :: & rTol + rTol = num%rTol_crystalliteState converged = all(abs(residuum) <= max(atol, rtol*abs(state))) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index af1108405..bd855e46f 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1325,18 +1325,19 @@ end function rhoDotFlux ! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! that sum up to a total of 1 are considered, all others are set to zero. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) +module subroutine plastic_nonlocal_updateCompatibility(orientation,ce) type(tRotationContainer), dimension(:), intent(in) :: & orientation ! crystal orientation integer, intent(in) :: & - ph, & - ip, & - el + ce integer :: & n, & ! neighbor index + ph, & en, & + ip, & + el, & neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor neighbor_me, & @@ -1344,17 +1345,21 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el) ns, & ! number of active slip systems s1, & ! slip system index (en) s2 ! slip system index (my neighbor) - real(pREAL), dimension(2,param(ph)%sum_N_sl,param(ph)%sum_N_sl,nIPneighbors) :: & + real(pREAL), dimension(2,param(material_ID_phase(1,ce))%sum_N_sl,param(material_ID_phase(1,ce))%sum_N_sl,nIPneighbors) :: & my_compatibility ! my_compatibility for current element and ip real(pREAL) :: & my_compatibilitySum, & thresholdValue, & nThresholdValues - logical, dimension(param(ph)%sum_N_sl) :: & + logical, dimension(param(material_ID_phase(1,ce))%sum_N_sl) :: & belowThreshold type(tRotation) :: mis + ph = material_ID_phase(1,ce) + el = (ce-1)/discretization_nIPs + 1 + ip = modulo(ce-1,discretization_nIPs) + 1 + associate(prm => param(ph)) ns = prm%sum_N_sl