diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index abbcce04a..240688a8c 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -5,7 +5,6 @@ !-------------------------------------------------------------------------------------------------- module CPFEM use prec - use FEsolving use math use rotations use YAML_types @@ -197,11 +196,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) else validCalculation - FEsolving_execElem = elCP - FEsolving_execIP = ip if (debugCPFEM%extensive) & print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip - call materialpoint_stressAndItsTangent(dt) + call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP]) terminalIllness: if (terminallyIll) then diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 44b93d1cb..5a500875d 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -6,7 +6,6 @@ module CPFEM2 use prec use config - use FEsolving use math use rotations use YAML_types diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index ea7430c6b..0ad68445c 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -176,7 +176,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & use DAMASK_interface use config use YAML_types - use FEsolving use discretization_marc use homogenization use CPFEM diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 deleted file mode 100644 index 3fc1482d3..000000000 --- a/src/FEsolving.f90 +++ /dev/null @@ -1,15 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief global variables for flow control -!-------------------------------------------------------------------------------------------------- -module FEsolving - - implicit none - public - - integer, dimension(2) :: & - FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element - FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP - -end module FEsolving diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 08e7b9c1c..d8ab6390d 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -13,7 +13,6 @@ #include "math.f90" #include "quaternions.f90" #include "rotations.f90" -#include "FEsolving.f90" #include "element.f90" #include "HDF5_utilities.f90" #include "results.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7e380f8cd..e65ce864d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -16,7 +16,6 @@ module constitutive use parallelization use HDF5_utilities use DAMASK_interface - use FEsolving use results implicit none @@ -65,10 +64,6 @@ module constitutive real(pReal), dimension(:,:,:,:,:), allocatable, public :: & crystallite_partitionedF !< def grad to be reached at end of homog inc - logical, dimension(:,:,:), allocatable :: & - crystallite_converged !< convergence flag - - type :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data end type @@ -186,10 +181,10 @@ module constitutive ! == cleaned:end =================================================================================== - module function crystallite_stress(dt,co,ip,el) + module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt integer, intent(in) :: co, ip, el - logical :: crystallite_stress + logical :: converged_ end function crystallite_stress module function constitutive_homogenizedC(co,ip,el) result(C) @@ -873,10 +868,8 @@ subroutine crystallite_init source = crystallite_partitionedF) allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal) - allocate(crystallite_orientation(cMax,iMax,eMax)) - allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) @@ -940,8 +933,8 @@ subroutine crystallite_init flush(IO_STDOUT) !$OMP PARALLEL DO PRIVATE(ph,me) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - do ip = FEsolving_execIP(1), FEsolving_execIP(2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + do el = 1, size(material_phaseMemberAt,3) + do ip = 1, size(material_phaseMemberAt,2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) @@ -967,14 +960,14 @@ subroutine crystallite_init crystallite_partitionedF0 = crystallite_F0 crystallite_partitionedF = crystallite_F0 - call crystallite_orientations() !$OMP PARALLEL DO PRIVATE(ph,me) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - do ip = FEsolving_execIP(1),FEsolving_execIP(2) + do el = 1, size(material_phaseMemberAt,3) + do ip = 1, size(material_phaseMemberAt,2) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) + call crystallite_orientations(co,ip,el) call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el) ! update dependent state variables to be consistent with basic states enddo enddo @@ -1089,7 +1082,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) el !< counter in element loop integer :: & o, & - p, pp, m + p, ph, me real(pReal), dimension(3,3) :: devNull, & invSubFp0,invSubFi0,invFp,invFi, & @@ -1109,19 +1102,19 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) real(pReal), dimension(9,9):: temp_99 logical :: error - pp = material_phaseAt(co,el) - m = material_phaseMemberAt(co,ip,el) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & crystallite_Fe(1:3,1:3,co,ip,el), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el) + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & crystallite_S (1:3,1:3,co,ip,el), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me), & co,ip,el) - invFp = math_inv33(constitutive_mech_Fp(pp)%data(1:3,1:3,m)) - invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m)) + invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) + invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el)) invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el)) @@ -1150,7 +1143,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & crystallite_S (1:3,1:3,co,ip,el), & - constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el) + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS !-------------------------------------------------------------------------------------------------- @@ -1210,34 +1203,20 @@ end function crystallite_stressTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- -subroutine crystallite_orientations +subroutine crystallite_orientations(co,ip,el) - integer & + integer, intent(in) :: & co, & !< counter in integration point component loop ip, & !< counter in integration point loop el !< counter in element loop - !$OMP PARALLEL DO - do el = FEsolving_execElem(1),FEsolving_execElem(2) - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el)))) - enddo; enddo; enddo - !$OMP END PARALLEL DO + call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el)))) + + if (plasticState(material_phaseAt(1,el))%nonlocal) & + call plastic_nonlocal_updateCompatibility(crystallite_orientation, & + phase_plasticityInstance(material_phaseAt(1,el)),ip,el) - nonlocalPresent: if (any(plasticState%nonlocal)) then - !$OMP PARALLEL DO - do el = FEsolving_execElem(1),FEsolving_execElem(2) - if (plasticState(material_phaseAt(1,el))%nonlocal) then - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - call plastic_nonlocal_updateCompatibility(crystallite_orientation, & - phase_plasticityInstance(material_phaseAt(1,el)),ip,el) - enddo - endif - enddo - !$OMP END PARALLEL DO - endif nonlocalPresent end subroutine crystallite_orientations @@ -1268,7 +1247,7 @@ end function crystallite_push33ToRef !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateSourceState(co,ip,el) +function integrateSourceState(co,ip,el) result(broken) integer, intent(in) :: & el, & !< element index in element loop @@ -1288,12 +1267,13 @@ subroutine integrateSourceState(co,ip,el) r ! state residuum real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & - broken + broken, converged_ ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) + converged_ = .true. broken = constitutive_thermal_collectDotState(ph,me) broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me) if(broken) return @@ -1328,19 +1308,20 @@ subroutine integrateSourceState(co,ip,el) - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * crystallite_subdt(co,ip,el) sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & - r(1:size_so(so)) - crystallite_converged(co,ip,el) = & - crystallite_converged(co,ip,el) .and. converged(r(1:size_so(so)), & - sourceState(ph)%p(so)%state(1:size_so(so),me), & - sourceState(ph)%p(so)%atol(1:size_so(so))) + converged_ = converged_ .and. converged(r(1:size_so(so)), & + sourceState(ph)%p(so)%state(1:size_so(so),me), & + sourceState(ph)%p(so)%atol(1:size_so(so))) enddo - if(crystallite_converged(co,ip,el)) then + if(converged_) then broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me) exit iteration endif enddo iteration + broken = broken .or. .not. converged_ + contains @@ -1364,7 +1345,7 @@ subroutine integrateSourceState(co,ip,el) end function damper -end subroutine integrateSourceState +end function integrateSourceState !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 7a2224ede..de6f2ae9f 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -951,7 +951,7 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el) +function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t @@ -1004,11 +1004,7 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el) - plasticState(ph)%dotState (1:size_pl,me) * Delta_t plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%state(1:size_pl,me) & - r(1:size_pl) - crystallite_converged(co,ip,el) = converged(r(1:size_pl), & - plasticState(ph)%state(1:size_pl,me), & - plasticState(ph)%atol(1:size_pl)) - - if(crystallite_converged(co,ip,el)) then + if (converged(r(1:size_pl),plasticState(ph)%state(1:size_pl,me),plasticState(ph)%atol(1:size_pl))) then broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) exit iteration @@ -1016,7 +1012,6 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el) enddo iteration - contains !-------------------------------------------------------------------------------------------------- @@ -1039,13 +1034,13 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el) end function damper -end subroutine integrateStateFPI +end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el) +function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t @@ -1075,15 +1070,14 @@ subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el) if(broken) return broken = integrateStress(F,Delta_t,co,ip,el) - crystallite_converged(co,ip,el) = .not. broken -end subroutine integrateStateEuler +end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) +function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t @@ -1123,24 +1117,22 @@ subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) if(broken) return + broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, & + plasticState(ph)%state(1:sizeDotState,me), & + plasticState(ph)%atol(1:sizeDotState)) - sizeDotState = plasticState(ph)%sizeDotState - crystallite_converged(co,ip,el) = converged(residuum_plastic(1:sizeDotState) & - + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, & - plasticState(ph)%state(1:sizeDotState,me), & - plasticState(ph)%atol(1:sizeDotState)) - -end subroutine integrateStateAdaptiveEuler +end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el) +function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t integer, intent(in) :: co,ip,el + logical :: broken real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1153,19 +1145,20 @@ subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el) real(pReal), dimension(4), parameter :: & B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] - call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C) + broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C) -end subroutine integrateStateRK4 +end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) +function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t integer, intent(in) :: co,ip,el + logical :: broken real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1185,16 +1178,16 @@ subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) + broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) -end subroutine integrateStateRKCK45 +end function integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) +function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken) real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in) :: Delta_t @@ -1205,15 +1198,14 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) el, & !< element index in element loop ip, & !< integration point index in ip loop co !< grain index in grain loop + logical :: broken - integer :: & + integer :: & stage, & ! stage index in integration stage loop n, & ph, & me, & sizeDotState - logical :: & - broken real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState @@ -1266,10 +1258,8 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) if(broken) return broken = integrateStress(F,Delta_t,co,ip,el) - crystallite_converged(co,ip,el) = .not. broken - -end subroutine integrateStateRK +end function integrateStateRK !-------------------------------------------------------------------------------------------------- @@ -1479,15 +1469,14 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function crystallite_stress(dt,co,ip,el) +module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt integer, intent(in) :: & co, & ip, & el - - logical :: crystallite_stress + logical :: converged_ real(pReal) :: & formerSubStep @@ -1519,7 +1508,7 @@ module function crystallite_stress(dt,co,ip,el) subFrac = 0.0_pReal subStep = 1.0_pReal/num%subStepSizeCryst todo = .true. - crystallite_converged(co,ip,el) = .false. ! pretend failed step of 1/subStepSizeCryst + converged_ = .false. ! pretend failed step of 1/subStepSizeCryst todo = .true. NiterationCrystallite = 0 @@ -1528,7 +1517,7 @@ module function crystallite_stress(dt,co,ip,el) !-------------------------------------------------------------------------------------------------- ! wind forward - if (crystallite_converged(co,ip,el)) then + if (converged_) then formerSubStep = subStep subFrac = subFrac + subStep subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep) @@ -1579,17 +1568,13 @@ module function crystallite_stress(dt,co,ip,el) math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) crystallite_subdt(co,ip,el) = subStep * dt - crystallite_converged(co,ip,el) = .false. - call integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),& - crystallite_subdt(co,ip,el),co,ip,el) - call integrateSourceState(co,ip,el) + converged_ = .not. integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),& + crystallite_subdt(co,ip,el),co,ip,el) + converged_ = converged_ .and. .not. integrateSourceState(co,ip,el) endif enddo cutbackLooping -! return whether converged or not - crystallite_stress = crystallite_converged(co,ip,el) - end function crystallite_stress end submodule constitutive_mech diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 1b3700c14..48ad5b7e1 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -19,7 +19,6 @@ module discretization_grid use results use discretization use geometry_plastic_nonlocal - use FEsolving implicit none private @@ -117,9 +116,6 @@ subroutine discretization_grid_init(restart) (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process worldrank+1==worldsize)) - FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements - FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP - !-------------------------------------------------------------------------------------------------- ! store geometry information for post processing if(.not. restart) then diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index cdf806b35..003f568c6 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -18,7 +18,6 @@ module grid_mech_FEM use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index ebaaf3b55..9bc36165f 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -18,7 +18,6 @@ module grid_mech_spectral_basic use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization_grid diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 9f2a17c97..7160c1adc 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -18,7 +18,6 @@ module grid_mech_spectral_polarisation use math use rotations use spectral_utilities - use FEsolving use config use homogenization use discretization_grid diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index c0c84233d..e8bae223a 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -810,9 +810,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& print'(/,a)', ' ... evaluating constitutive response ......................................' flush(IO_STDOUT) - homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field - call materialpoint_stressAndItsTangent(timeinc) ! calculate P field + call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 27fdb6064..ebf5fd50d 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -11,7 +11,6 @@ module homogenization use math use material use constitutive - use FEsolving use discretization use thermal_isothermal use thermal_conduction @@ -144,27 +143,29 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- !> @brief parallelized calculation of stress and corresponding tangent at material points !-------------------------------------------------------------------------------------------------- -subroutine materialpoint_stressAndItsTangent(dt) +subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem) real(pReal), intent(in) :: dt !< time increment + integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP integer :: & NiterationHomog, & NiterationMPstate, & ip, & !< integration point number el, & !< element number - myNgrains, co, ce + myNgrains, co, ce, ho real(pReal) :: & subFrac, & subStep logical :: & - requested, & converged logical, dimension(2) :: & doneAndHappy -!$OMP PARALLEL DO PRIVATE(ce,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,requested,doneAndHappy) +!$OMP PARALLEL DO PRIVATE(ce,ho,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,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) !-------------------------------------------------------------------------------------------------- @@ -174,21 +175,19 @@ subroutine materialpoint_stressAndItsTangent(dt) subFrac = 0.0_pReal converged = .false. ! pretend failed step ... subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation - requested = .true. ! everybody requires calculation - if (homogState(material_homogenizationAt(el))%sizeState > 0) & - homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - homogState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el)) + if (homogState(ho)%sizeState > 0) & + homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + homogState(ho)%State0( :,material_homogenizationMemberAt(ip,el)) + + if (damageState(ho)%sizeState > 0) & + damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + damageState(ho)%State0( :,material_homogenizationMemberAt(ip,el)) - if (damageState(material_homogenizationAt(el))%sizeState > 0) & - damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - damageState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el)) NiterationHomog = 0 cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) - myNgrains = homogenization_Nconstituents(material_homogenizationAt(el)) - if (converged) then subFrac = subFrac + subStep subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration @@ -198,22 +197,20 @@ subroutine materialpoint_stressAndItsTangent(dt) ! wind forward grain starting point call constitutive_windForward(ip,el) - if(homogState(material_homogenizationAt(el))%sizeState > 0) & - homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - homogState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el)) - if(damageState(material_homogenizationAt(el))%sizeState > 0) & - damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = & - damageState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el)) + if(homogState(ho)%sizeState > 0) & + homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + homogState(ho)%State (:,material_homogenizationMemberAt(ip,el)) + if(damageState(ho)%sizeState > 0) & + damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = & + damageState(ho)%State (:,material_homogenizationMemberAt(ip,el)) endif steppingNeeded - else if ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep ! cutback makes no sense - if (.not. terminallyIll) then ! so first signals terminally ill... + if (.not. terminallyIll) & ! so first signals terminally ill... print*, ' Integration point ', ip,' at element ', el, ' terminally ill' - endif terminallyIll = .true. ! ...and kills all others else ! cutback makes sense subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback @@ -221,23 +218,19 @@ subroutine materialpoint_stressAndItsTangent(dt) call crystallite_restore(ip,el,subStep < 1.0_pReal) call constitutive_restore(ip,el) - if(homogState(material_homogenizationAt(el))%sizeState > 0) & - homogState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = & - homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) - if(damageState(material_homogenizationAt(el))%sizeState > 0) & - damageState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = & - damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) + if(homogState(ho)%sizeState > 0) & + homogState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = & + homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) + if(damageState(ho)%sizeState > 0) & + damageState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = & + damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) endif endif - if (subStep > num%subStepMinHomog) then - requested = .true. - doneAndHappy = [.false.,.true.] - endif - + if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] NiterationMPstate = 0 - convergenceLooping: do while (.not. terminallyIll .and. requested & + convergenceLooping: do while (.not. terminallyIll & .and. .not. doneAndHappy(1) & .and. NiterationMPstate < num%nMPstate) NiterationMPstate = NiterationMPstate + 1 @@ -245,7 +238,7 @@ subroutine materialpoint_stressAndItsTangent(dt) !-------------------------------------------------------------------------------------------------- ! deformation partitioning - if(requested .and. .not. doneAndHappy(1)) then ! requested but not yet done + if (.not. doneAndHappy(1)) then ce = (el-1)*discretization_nIPs + ip call mech_partition(homogenization_F0(1:3,1:3,ce) & + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))& @@ -255,10 +248,7 @@ subroutine materialpoint_stressAndItsTangent(dt) do co = 1, myNgrains converged = converged .and. crystallite_stress(dt*subStep,co,ip,el) enddo - endif - - if (requested .and. .not. doneAndHappy(1)) then if (.not. converged) then doneAndHappy = [.true.,.false.] else @@ -281,10 +271,14 @@ subroutine materialpoint_stressAndItsTangent(dt) !$OMP END PARALLEL DO if (.not. terminallyIll ) then - call crystallite_orientations() ! calculate crystal orientations - !$OMP PARALLEL DO + !$OMP PARALLEL DO PRIVATE(ho,myNgrains) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) + ho = material_homogenizationAt(el) + myNgrains = homogenization_Nconstituents(ho) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) + do co = 1, myNgrains + call crystallite_orientations(co,ip,el) + enddo call mech_homogenize(ip,el) enddo IpLooping3 enddo elementLooping3 diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 56f1e554f..e4499e9b7 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -128,35 +128,35 @@ module subroutine mech_homogenize(ip,el) integer, intent(in) :: & ip, & !< integration point el !< element number - integer :: c,m + integer :: co,ce real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) - m = (el-1)* discretization_nIPs + ip + ce = (el-1)* discretization_nIPs + ip chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,m) = crystallite_P(1:3,1:3,1,ip,el) - homogenization_dPdF(1:3,1:3,1:3,1:3,m) = crystallite_stressTangent(1,ip,el) + homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,m), & - homogenization_dPdF(1:3,1:3,1:3,1:3,m),& + homogenization_P(1:3,1:3,ce), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& - homogenization_P(1:3,1:3,m), & - homogenization_dPdF(1:3,1:3,1:3,1:3,m),& + homogenization_P(1:3,1:3,ce), & + homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index ca0b54b73..675e57bd3 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -12,7 +12,6 @@ module discretization_marc use DAMASK_interface use IO use config - use FEsolving use element use discretization use geometry_plastic_nonlocal @@ -89,9 +88,6 @@ subroutine discretization_marc_init if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element') if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [1,nElems] - FEsolving_execIP = [1,elem%nIPs] - allocate(cellNodeDefinition(elem%nNodes-1)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 1e353892c..7369520c1 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -15,7 +15,6 @@ program DAMASK_mesh use IO use math use CPFEM2 - use FEsolving use config use discretization_mesh use FEM_Utilities diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index cb81f1f0c..2f3633e11 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -160,7 +160,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) print'(/,a)', ' ... evaluating constitutive response ......................................' - call materialpoint_stressAndItsTangent(timeinc) ! calculate P field + call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field cutBack = .false. ! reset cutBack status diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 7dbd05e46..21c5feace 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -18,7 +18,6 @@ module discretization_mesh use config use discretization use results - use FEsolving use FEM_quadrature use YAML_types use prec @@ -30,7 +29,7 @@ module discretization_mesh mesh_Nboundaries, & mesh_NcpElemsGlobal - integer :: & + integer, public, protected :: & mesh_NcpElems !< total number of CP elements in mesh !!!! BEGIN DEPRECATED !!!!! @@ -174,9 +173,6 @@ subroutine discretization_mesh_init(restart) if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element') if (debug_ip < 1 .or. debug_ip > mesh_maxNips) call IO_error(602,ext_msg='IP') - FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements - FEsolving_execIP = [1,mesh_maxNips] - allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) call discretization_init(materialAt,&