Merge branch 'solver-cleanup' into 'development'

Solver cleanup

See merge request damask/DAMASK!888

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}")

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
 

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]
 

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]
 

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

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
 

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

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)
+  !$OMP PARALLEL DO PRIVATE(ho)
-  elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
+  do ce = cell_start, cell_end
-    IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
-      ce = (el-1)*discretization_nIPs + ip
       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
-  end do elementLooping3
   !$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
 

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

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)
+    module subroutine plastic_nonlocal_updateCompatibility(orientation,ce)
-      integer, intent(in) :: &
+      integer, intent(in) :: ce
-        ph, &
-        ip, &
-        el
         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)
+  !$OMP PARALLEL DO PRIVATE(ph,en)
-  do el = 1, discretization_Nelems
+  do ce = 1, size(material_ID_homogenization)
-    do ip = 1, discretization_nIPs
+    do co = 1,homogenization_Nconstituents(material_ID_homogenization(ce))
-      ce = (el-1)*discretization_nIPs + ip
+      ph = material_ID_phase(co,ce)
-      do co = 1,homogenization_Nconstituents(material_ID_homogenization(ce))
+      en = material_entry_phase(co,ce)
-        en = material_entry_phase(co,ce)
+      call crystallite_orientations(co,ce)
-        ph = material_ID_phase(co,ce)
+      call plastic_dependentState(ph,en)                                                          ! update dependent state variables to be consistent with basic states
-        call crystallite_orientations(co,ip,el)
-        call plastic_dependentState(ph,en)                                                          ! update dependent state variables to be consistent with basic states
-     end do
     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
+    co, &
-    ip, &                                                                                           !< counter in integration point loop
+    ce
-    el                                                                                              !< counter in element loop
 
   integer :: ph, en
 
 
-  ph = material_ID_phase(co,(el-1)*discretization_nIPs + ip)
+  ph = material_ID_phase(co,ce)
-  en = material_entry_phase(co,(el-1)*discretization_nIPs + ip)
+  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) &
+  if (plasticState(material_ID_phase(1,ce))%nonlocal) call plastic_nonlocal_updateCompatibility(phase_O,ce)
-    call plastic_nonlocal_updateCompatibility(phase_O,material_ID_phase(1,(el-1)*discretization_nIPs + ip),ip,el)
 
 
 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)))

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, &
+    ce
-    ip, &
-    el
 
   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