diff --git a/examples/config/phase/damage/isobrittle_generic.yaml b/examples/config/phase/damage/isobrittle_generic.yaml index 98c9b9763..a7c30ff1e 100644 --- a/examples/config/phase/damage/isobrittle_generic.yaml +++ b/examples/config/phase/damage/isobrittle_generic.yaml @@ -1,7 +1,5 @@ type: isobrittle W_crit: 1400000.0 -m: 1.0 -isoBrittle_atol: 0.01 output: [f_phi] diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 213b56a54..aa532859a 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -189,9 +189,11 @@ 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 - if (debugCPFEM%extensive) & - print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip - call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP]) + if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip + call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP]) + if (.not. terminallyIll) & + call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP]) + terminalIllness: if (terminallyIll) then diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 77678137d..111e9fece 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -261,7 +261,7 @@ subroutine grid_mechanical_FEM_init F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F @@ -490,8 +490,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i BCTol err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ - divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol) - BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol) + divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol) + BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pReal)) & .or. terminallyIll) then @@ -502,8 +502,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i reason = 0 endif -!-------------------------------------------------------------------------------------------------- -! report print'(1/,a)', ' ... reporting .............................................................' print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 83b961023..c4fce61a7 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -211,7 +211,7 @@ subroutine grid_mechanical_spectral_basic_init F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -429,8 +429,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm divTol, & BCTol - divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol) - BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol) + divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol) + BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pReal)) & .or. terminallyIll) then @@ -441,8 +441,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm reason = 0 endif -!-------------------------------------------------------------------------------------------------- -! report print'(1/,a)', ' ... reporting .............................................................' print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 07edc84b5..a52554e3c 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -237,7 +237,7 @@ subroutine grid_mechanical_spectral_polarisation_init F_tau_lastInc = 2.0_pReal*F_lastInc endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -487,9 +487,9 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm divTol, & BCTol - curlTol = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol ,num%eps_curl_atol) - divTol = max(maxval(abs(P_av)) *num%eps_div_rtol ,num%eps_div_atol) - BCTol = max(maxval(abs(P_av)) *num%eps_stress_rtol,num%eps_stress_atol) + curlTol = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol, num%eps_curl_atol) + divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol) + BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol) if ((totalIter >= num%itmin .and. all([err_div/divTol, err_curl/curlTol, err_BC/BCTol] < 1.0_pReal)) & .or. terminallyIll) then @@ -500,8 +500,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm reason = 0 endif -!-------------------------------------------------------------------------------------------------- -! report print'(1/,a)', ' ... reporting .............................................................' print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' @@ -542,13 +540,13 @@ subroutine formResidual(in, FandF_tau, & !--------------------------------------------------------------------------------------------------- F => FandF_tau(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) + XG_RANGE,YG_RANGE,ZG_RANGE) F_tau => FandF_tau(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) + XG_RANGE,YG_RANGE,ZG_RANGE) residual_F => residuum(1:3,1:3,1,& - X_RANGE, Y_RANGE, Z_RANGE) + X_RANGE, Y_RANGE, Z_RANGE) residual_F_tau => residuum(1:3,1:3,2,& - X_RANGE, Y_RANGE, Z_RANGE) + X_RANGE, Y_RANGE, Z_RANGE) F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 061cc9d74..bd37e2654 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -815,10 +815,14 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field - call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field + call homogenization_mechanical_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field + if (.not. terminallyIll) & + call homogenization_thermal_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + if (.not. terminallyIll) & + call homogenization_mechanical_response2(timeinc,[1,1],[1,product(grid(1:2))*grid3]) 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 + P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & ' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal diff --git a/src/homogenization.f90 b/src/homogenization.f90 index f34834272..fc5dd7bda 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -101,8 +101,8 @@ module homogenization integer, intent(in) :: ce end subroutine damage_partition - module subroutine mechanical_homogenize(dt,ce) - real(pReal), intent(in) :: dt + module subroutine mechanical_homogenize(Delta_t,ce) + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & ce !< cell end subroutine mechanical_homogenize @@ -178,7 +178,9 @@ module homogenization public :: & homogenization_init, & - materialpoint_stressAndItsTangent, & + homogenization_mechanical_response, & + homogenization_mechanical_response2, & + homogenization_thermal_response, & homogenization_mu_T, & homogenization_K_T, & homogenization_f_T, & @@ -227,24 +229,24 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- -!> @brief parallelized calculation of stress and corresponding tangent at material points +!> @brief !-------------------------------------------------------------------------------------------------- -subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem) - real(pReal), intent(in) :: dt !< time increment + real(pReal), intent(in) :: Delta_t !< time increment integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP integer :: & NiterationMPstate, & ip, & !< integration point number el, & !< element number - co, ce, ho, en, ph + co, ce, ho, en logical :: & converged logical, dimension(2) :: & doneAndHappy - !$OMP PARALLEL - !$OMP DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy) + + !$OMP PARALLEL DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) do ip = FEsolving_execIP(1),FEsolving_execIP(2) @@ -266,65 +268,89 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE NiterationMPstate = NiterationMPstate + 1 call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) - converged = .true. - do co = 1, homogenization_Nconstituents(ho) - converged = converged .and. crystallite_stress(dt,co,ip,el) - enddo - + converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) if (converged) then - doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ce) + doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) converged = all(doneAndHappy) else doneAndHappy = [.true.,.false.] endif - enddo convergenceLooping + + converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) + if (.not. converged) then - if (.not. terminallyIll) print*, ' Integration point ', ip,' at element ', el, ' terminally ill' + if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' terminallyIll = .true. endif enddo enddo - !$OMP END DO + !$OMP END PARALLEL DO - if (.not. terminallyIll) then - !$OMP DO PRIVATE(ho,ph,ce) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - if (terminallyIll) continue - 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_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' - terminallyIll = .true. ! ...and kills all others - endif - enddo +end subroutine homogenization_mechanical_response + + +!-------------------------------------------------------------------------------------------------- +!> @brief +!-------------------------------------------------------------------------------------------------- +subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem) + + real(pReal), intent(in) :: Delta_t !< time increment + integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + integer :: & + ip, & !< integration point number + el, & !< element number + co, ce, ho + + + !$OMP PARALLEL DO PRIVATE(ho,ce) + do el = FEsolving_execElem(1),FEsolving_execElem(2) + if (terminallyIll) continue + 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) + if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then + if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' + terminallyIll = .true. + endif enddo enddo - !$OMP END DO + enddo + !$OMP END PARALLEL DO - !$OMP 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 - ho = material_homogenizationID(ce) - do co = 1, homogenization_Nconstituents(ho) - call crystallite_orientations(co,ip,el) - enddo - call mechanical_homogenize(dt,ce) - enddo IpLooping3 - enddo elementLooping3 - !$OMP END DO - else - print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' - endif - !$OMP END PARALLEL +end subroutine homogenization_thermal_response -end subroutine materialpoint_stressAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief +!-------------------------------------------------------------------------------------------------- +subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolving_execElem) + + real(pReal), intent(in) :: Delta_t !< time increment + integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP + 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 + ho = material_homogenizationID(ce) + do co = 1, homogenization_Nconstituents(ho) + call crystallite_orientations(co,ip,el) + enddo + call mechanical_homogenize(Delta_t,ce) + enddo IpLooping3 + enddo elementLooping3 + !$OMP END PARALLEL DO + + +end subroutine homogenization_mechanical_response2 !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index 7ea1f74e9..0138d2468 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -123,21 +123,21 @@ end subroutine mechanical_partition !-------------------------------------------------------------------------------------------------- !> @brief Average P and dPdF from the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_homogenize(dt,ce) +module subroutine mechanical_homogenize(Delta_t,ce) - real(pReal), intent(in) :: dt + real(pReal), intent(in) :: Delta_t integer, intent(in) :: ce integer :: co homogenization_P(1:3,1:3,ce) = phase_P(1,ce) - homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) & + phase_P(co,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) & - + phase_mechanical_dPdF(dt,co,ce) + + phase_mechanical_dPdF(Delta_t,co,ce) enddo homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) & diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 6765d3d0d..5bb6d5f44 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -27,17 +27,17 @@ module FEM_utilities logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill integer, public, parameter :: maxFields = 6 integer, public :: nActiveFields = 0 - + !-------------------------------------------------------------------------------------------------- ! grid related information information real(pReal), public :: wgt !< weighting factor 1/Nelems - + !-------------------------------------------------------------------------------------------------- ! field labels information character(len=*), parameter, public :: & FIELD_MECH_label = 'mechanical' - + enum, bind(c); enumerator :: & FIELD_UNDEFINED_ID, & FIELD_MECH_ID @@ -48,7 +48,7 @@ module FEM_utilities COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Z_ID end enum - + !-------------------------------------------------------------------------------------------------- ! variables controlling debugging logical :: & @@ -57,23 +57,23 @@ module FEM_utilities !-------------------------------------------------------------------------------------------------- ! derived types type, public :: tSolutionState !< return type of solution from FEM solver variants - logical :: converged = .true. - logical :: stagConverged = .true. + logical :: converged = .true. + logical :: stagConverged = .true. integer :: iterationsNeeded = 0 end type tSolutionState - + type, public :: tComponentBC integer(kind(COMPONENT_UNDEFINED_ID)) :: ID real(pReal), allocatable, dimension(:) :: Value - logical, allocatable, dimension(:) :: Mask + logical, allocatable, dimension(:) :: Mask end type tComponentBC - + type, public :: tFieldBC integer(kind(FIELD_UNDEFINED_ID)) :: ID integer :: nComponents = 0 type(tComponentBC), allocatable :: componentBC(:) end type tFieldBC - + type, public :: tLoadCase real(pReal) :: time = 0.0_pReal !< length of increment integer :: incs = 0, & !< number of increments @@ -83,7 +83,7 @@ module FEM_utilities integer, allocatable, dimension(:) :: faceID type(tFieldBC), allocatable, dimension(:) :: fieldBC end type tLoadCase - + public :: & FEM_utilities_init, & utilities_constitutiveResponse, & @@ -94,14 +94,14 @@ module FEM_utilities COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Z_ID -contains +contains !ToDo: use functions in variable call !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, sets debug flags !-------------------------------------------------------------------------------------------------- subroutine FEM_utilities_init - + character(len=pStringLen) :: petsc_optionsOrder class(tNode), pointer :: & num_mesh, & @@ -113,7 +113,7 @@ subroutine FEM_utilities_init PetscErrorCode :: ierr print'(/,a)', ' <<<+- FEM_utilities init -+>>>' - + num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2) @@ -141,7 +141,7 @@ subroutine FEM_utilities_init write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr) CHKERRQ(ierr) - + wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal) @@ -152,20 +152,21 @@ end subroutine FEM_utilities_init !> @brief calculates constitutive response !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) - + real(pReal), intent(in) :: timeinc !< loading time logical, intent(in) :: forwardData !< age results - + real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress - + PetscErrorCode :: ierr print'(/,a)', ' ... evaluating constitutive response ......................................' - call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field + call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field + if (.not. terminallyIll) & + call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) + cutBack = .false. - cutBack = .false. ! reset cutBack status - P_av = sum(homogenization_P,dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) @@ -198,12 +199,12 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa do dof = offset+comp+1, offset+numDof, numComp localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc enddo - enddo + enddo call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr) call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr) call VecAssemblyEnd (localVec, ierr); CHKERRQ(ierr) if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr) - + end subroutine utilities_projectBCValues end module FEM_utilities diff --git a/src/phase.f90 b/src/phase.f90 index e396e5661..d90f41ccc 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -108,9 +108,13 @@ module phase logical, intent(in) :: includeL end subroutine mechanical_restore + module subroutine damage_restore(ce) + integer, intent(in) :: ce + end subroutine damage_restore - module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) - real(pReal), intent(in) :: dt + + module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF) + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & !< counter in constituent loop ce @@ -164,8 +168,8 @@ module phase real(pReal) :: dot_T end function thermal_dot_T - module function damage_phi(ph,me) result(phi) - integer, intent(in) :: ph,me + module function damage_phi(ph,en) result(phi) + integer, intent(in) :: ph,en real(pReal) :: phi end function damage_phi @@ -209,29 +213,27 @@ module phase ! == cleaned:end =================================================================================== - module function thermal_stress(Delta_t,ph,en) result(converged_) + module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) real(pReal), intent(in) :: Delta_t integer, intent(in) :: ph, en logical :: converged_ - end function thermal_stress + end function phase_thermal_constitutive - module function integrateDamageState(dt,co,ce) result(broken) - real(pReal), intent(in) :: dt - integer, intent(in) :: & - ce, & - co - logical :: broken - end function integrateDamageState - - module function crystallite_stress(dt,co,ip,el) result(converged_) - real(pReal), intent(in) :: dt + module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) + real(pReal), intent(in) :: Delta_t integer, intent(in) :: co, ip, el logical :: converged_ - end function crystallite_stress + end function phase_damage_constitutive - !ToDo: Try to merge the all stiffness functions + module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_) + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: co, ip, el + logical :: converged_ + end function phase_mechanical_constitutive + + !ToDo: Merge all the stiffness functions module function phase_homogenizedC(ph,en) result(C) integer, intent(in) :: ph, en real(pReal), dimension(6,6) :: C @@ -271,8 +273,8 @@ module phase end subroutine plastic_dependentState - module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) - integer, intent(in) :: ph, me + module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en) + integer, intent(in) :: ph, en real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -315,14 +317,14 @@ module phase plastic_nonlocal_updateCompatibility, & converged, & crystallite_init, & - crystallite_stress, & - thermal_stress, & + phase_mechanical_constitutive, & + phase_thermal_constitutive, & + phase_damage_constitutive, & phase_mechanical_dPdF, & crystallite_orientations, & crystallite_push33ToRef, & phase_restartWrite, & phase_restartRead, & - integrateDamageState, & phase_thermal_setField, & phase_set_phi, & phase_P, & @@ -435,17 +437,9 @@ subroutine phase_restore(ce,includeL) logical, intent(in) :: includeL integer, intent(in) :: ce - integer :: & - co - - - do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) - if (damageState(material_phaseID(co,ce))%sizeState > 0) & - damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = & - damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce)) - enddo call mechanical_restore(ce,includeL) + call damage_restore(ce) end subroutine phase_restore @@ -545,10 +539,6 @@ subroutine crystallite_init() phases => config_material%get('phase') - do ph = 1, phases%length - if (damageState(ph)%sizeState > 0) allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack - enddo - print'(a42,1x,i10)', ' # of elements: ', eMax print'(a42,1x,i10)', ' # of integration points/element: ', iMax print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index d8ac5046a..172f6f22a 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -39,8 +39,8 @@ submodule(phase) damage end function isobrittle_init - module subroutine isobrittle_deltaState(C, Fe, ph, me) - integer, intent(in) :: ph,me + module subroutine isobrittle_deltaState(C, Fe, ph, en) + integer, intent(in) :: ph,en real(pReal), intent(in), dimension(3,3) :: & Fe real(pReal), intent(in), dimension(6,6) :: & @@ -48,8 +48,8 @@ submodule(phase) damage end subroutine isobrittle_deltaState - module subroutine anisobrittle_dotState(S, ph, me) - integer, intent(in) :: ph,me + module subroutine anisobrittle_dotState(S, ph, en) + integer, intent(in) :: ph,en real(pReal), intent(in), dimension(3,3) :: & S end subroutine anisobrittle_dotState @@ -125,6 +125,29 @@ module subroutine damage_init end subroutine damage_init +!-------------------------------------------------------------------------------------------------- +!> @brief calculate stress (P) +!-------------------------------------------------------------------------------------------------- +module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) + + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: & + co, & + ip, & + el + logical :: converged_ + + integer :: & + ph, en + + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) + en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + + converged_ = .not. integrateDamageState(Delta_t,ph,en) + +end function phase_damage_constitutive + + !-------------------------------------------------------------------------------------------------- !> @brief returns the degraded/modified elasticity matrix !-------------------------------------------------------------------------------------------------- @@ -144,6 +167,26 @@ module function phase_damage_C(C_homogenized,ph,en) result(C) end function phase_damage_C +!-------------------------------------------------------------------------------------------------- +!> @brief Restore data after homog cutback. +!-------------------------------------------------------------------------------------------------- +module subroutine damage_restore(ce) + + integer, intent(in) :: ce + + integer :: & + co + + + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + if (damageState(material_phaseID(co,ce))%sizeState > 0) & + damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = & + damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce)) + enddo + +end subroutine damage_restore + + !---------------------------------------------------------------------------------------------- !< @brief returns local part of nonlocal damage driving force !---------------------------------------------------------------------------------------------- @@ -173,23 +216,20 @@ module function phase_f_phi(phi,co,ce) result(f) 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,ce) result(broken) +function integrateDamageState(Delta_t,ph,en) result(broken) - real(pReal), intent(in) :: dt + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - ce, & - co + ph, & + en logical :: broken integer :: & NiterationState, & !< number of iterations in state loop - ph, & - me, & size_so real(pReal) :: & zeta @@ -199,8 +239,6 @@ module function integrateDamageState(dt,co,ce) result(broken) logical :: & converged_ - ph = material_phaseID(co,ce) - me = material_phaseEntry(co,ce) if (damageState(ph)%sizeState == 0) then broken = .false. @@ -208,37 +246,37 @@ module function integrateDamageState(dt,co,ce) result(broken) endif converged_ = .true. - broken = phase_damage_collectDotState(ph,me) + broken = phase_damage_collectDotState(ph,en) if(broken) return - size_so = damageState(ph)%sizeDotState - damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) & - + damageState(ph)%dotState (1:size_so,me) * dt - source_dotState(1:size_so,2) = 0.0_pReal + size_so = damageState(ph)%sizeDotState + damageState(ph)%state(1:size_so,en) = damageState(ph)%state0 (1:size_so,en) & + + damageState(ph)%dotState(1:size_so,en) * Delta_t + source_dotState(1:size_so,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) - source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me) + source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en) - broken = phase_damage_collectDotState(ph,me) + broken = phase_damage_collectDotState(ph,en) if(broken) exit iteration - zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) - damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta & + zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) + damageState(ph)%dotState(:,en) = damageState(ph)%dotState(:,en) * zeta & + source_dotState(1:size_so,1)* (1.0_pReal - zeta) - r(1:size_so) = damageState(ph)%state (1:size_so,me) & - - damageState(ph)%subState0(1:size_so,me) & - - damageState(ph)%dotState (1:size_so,me) * dt - damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so) + r(1:size_so) = damageState(ph)%state (1:size_so,en) & + - damageState(ph)%State0 (1:size_so,en) & + - damageState(ph)%dotState(1:size_so,en) * Delta_t + damageState(ph)%state(1:size_so,en) = damageState(ph)%state(1:size_so,en) - r(1:size_so) converged_ = converged_ .and. converged(r(1:size_so), & - damageState(ph)%state(1:size_so,me), & + damageState(ph)%state(1:size_so,en), & damageState(ph)%atol(1:size_so)) if(converged_) then - broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me) + broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en) exit iteration endif @@ -300,11 +338,11 @@ end subroutine damage_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function phase_damage_collectDotState(ph,me) result(broken) +function phase_damage_collectDotState(ph,en) result(broken) integer, intent(in) :: & ph, & - me !< counter in source loop + en !< counter in source loop logical :: broken @@ -315,11 +353,11 @@ function phase_damage_collectDotState(ph,me) result(broken) sourceType: select case (phase_damage(ph)) case (DAMAGE_ANISOBRITTLE_ID) sourceType - call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress? + call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! correct stress? end select sourceType - broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me))) + broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en))) endif @@ -358,11 +396,11 @@ end function phase_K_phi !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -function phase_damage_deltaState(Fe, ph, me) result(broken) +function phase_damage_deltaState(Fe, ph, en) result(broken) integer, intent(in) :: & ph, & - me + en real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer :: & @@ -379,13 +417,13 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) sourceType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) - broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me))) + call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en) + broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) if(.not. broken) then myOffset = damageState(ph)%offsetDeltaState mySize = damageState(ph)%sizeDeltaState - damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = & - damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me) + damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = & + damageState(ph)%state(myOffset + 1: myOffset + mySize,en) + damageState(ph)%deltaState(1:mySize,en) endif end select sourceType @@ -436,13 +474,13 @@ module subroutine phase_set_phi(phi,co,ce) end subroutine phase_set_phi -module function damage_phi(ph,me) result(phi) +module function damage_phi(ph,en) result(phi) - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en real(pReal) :: phi - phi = current(ph)%phi(me) + phi = current(ph)%phi(en) end function damage_phi diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 9984e5514..213d04f28 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -40,7 +40,7 @@ module function anisobrittle_init() result(mySources) phase, & sources, & src - integer :: Nmembers,p + integer :: Nmembers,ph integer, dimension(:), allocatable :: N_cl character(len=pStringLen) :: extmsg = '' @@ -56,12 +56,12 @@ module function anisobrittle_init() result(mySources) allocate(param(phases%length)) - do p = 1, phases%length - if(mySources(p)) then - phase => phases%get(p) - sources => phase%get('damage') + do ph = 1, phases%length + if(mySources(ph)) then + phase => phases%get(ph) + sources => phase%get('damage') - associate(prm => param(p)) + associate(prm => param(ph)) src => sources%get(1) N_cl = src%get_as1dInt('N_cl',defaultVal=emptyIntArray) @@ -92,17 +92,15 @@ 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_phaseID==p) - call phase_allocateState(damageState(p),Nmembers,1,1,0) - damageState(p)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) - if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' + Nmembers = count(material_phaseID==ph) + call phase_allocateState(damageState(ph),Nmembers,1,1,0) + damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) + if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' - end associate + end associate -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)') - endif + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)') + endif enddo @@ -112,10 +110,10 @@ end function anisobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -module subroutine anisobrittle_dotState(S, ph,me) +module subroutine anisobrittle_dotState(S, ph,en) integer, intent(in) :: & - ph,me + ph,en real(pReal), intent(in), dimension(3,3) :: & S @@ -126,15 +124,15 @@ module subroutine anisobrittle_dotState(S, ph,me) associate(prm => param(ph)) - damageState(ph)%dotState(1,me) = 0.0_pReal + damageState(ph)%dotState(1,en) = 0.0_pReal do i = 1, prm%sum_N_cl traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal + traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal - damageState(ph)%dotState(1,me) = damageState(ph)%dotState(1,me) & + damageState(ph)%dotState(1,en) = damageState(ph)%dotState(1,en) & + prm%dot_o / prm%s_crit(i) & * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & @@ -171,10 +169,10 @@ end subroutine anisobrittle_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) +module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en) integer, intent(in) :: & - ph,me + ph,en real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -193,7 +191,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) dLd_dTstar = 0.0_pReal associate(prm => param(ph)) do i = 1,prm%sum_N_cl - traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal + traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) if (abs(traction_d) > traction_crit + tol_math_check) then diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index aed6a91e9..066cab47e 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -13,7 +13,15 @@ submodule(phase:damage) isobrittle output end type tParameters - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) + type :: tIsobrittleState + real(pReal), pointer, dimension(:) :: & !< vectors along Nmembers + r_W !< ratio between actual and critical strain energy density + end type tIsobrittleState + + type(tParameters), allocatable, dimension(:) :: param !< containers of constitutive parameters (len Ninstances) + type(tIsobrittleState), allocatable, dimension(:) :: & + deltaState, & + state contains @@ -44,13 +52,15 @@ module function isobrittle_init() result(mySources) phases => config_material%get('phase') allocate(param(phases%length)) + allocate(state(phases%length)) + allocate(deltaState(phases%length)) do ph = 1, phases%length if(mySources(ph)) then - phase => phases%get(ph) - sources => phase%get('damage') + phase => phases%get(ph) + sources => phase%get('damage') - associate(prm => param(ph)) + associate(prm => param(ph), dlt => deltaState(ph), stt => state(ph)) src => sources%get(1) prm%W_crit = src%get_asFloat('W_crit') @@ -69,12 +79,14 @@ module function isobrittle_init() result(mySources) damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' - end associate + stt%r_W => damageState(ph)%state(1,:) + dlt%r_W => damageState(ph)%deltaState(1,:) -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)') - endif + end associate + + + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)') + endif enddo @@ -85,29 +97,27 @@ end function isobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine isobrittle_deltaState(C, Fe, ph,me) +module subroutine isobrittle_deltaState(C, Fe, ph,en) - integer, intent(in) :: ph,me + integer, intent(in) :: ph,en real(pReal), intent(in), dimension(3,3) :: & Fe real(pReal), intent(in), dimension(6,6) :: & C real(pReal), dimension(6) :: & - strain + epsilon real(pReal) :: & - strainenergy + r_W - strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) + epsilon = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) - associate(prm => param(ph)) - strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit - ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit + associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) + + r_W = 2.0_pReal*dot_product(epsilon,matmul(C,epsilon))/prm%W_crit + dlt%r_W(en) = merge(r_W - stt%r_W(en), 0.0_pReal, r_W > stt%r_W(en)) - damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), & - damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me), & - strainenergy > damageState(ph)%subState0(1,me)) end associate end subroutine isobrittle_deltaState @@ -124,14 +134,15 @@ module subroutine isobrittle_results(phase,group) integer :: o - associate(prm => param(phase), & - stt => damageState(phase)%state) + associate(prm => param(phase), stt => damageState(phase)%state) ! point to state and output r_W (is scalar, not 1D vector) + outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') - call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') + call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') ! Wrong, this is dimensionless end select enddo outputsLoop + end associate end subroutine isobrittle_results diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 7562c055f..8b04a0009 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -44,9 +44,9 @@ submodule(phase) mechanical interface - module subroutine eigendeformation_init(phases) + module subroutine eigen_init(phases) class(tNode), pointer :: phases - end subroutine eigendeformation_init + end subroutine eigen_init module subroutine elastic_init(phases) class(tNode), pointer :: phases @@ -197,8 +197,7 @@ module subroutine mechanical_init(materials,phases) Nmembers class(tNode), pointer :: & num_crystallite, & - material, & - constituents, & + phase, & mech @@ -303,7 +302,7 @@ module subroutine mechanical_init(materials,phases) end select - call eigendeformation_init(phases) + call eigen_init(phases) end subroutine mechanical_init @@ -977,9 +976,9 @@ end subroutine mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function crystallite_stress(dt,co,ip,el) result(converged_) +module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_) - real(pReal), intent(in) :: dt + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & ip, & @@ -1009,17 +1008,13 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) subState0 = plasticState(ph)%State0(:,en) - - if (damageState(ph)%sizeState > 0) & - damageState(ph)%subState0(:,en) = damageState(ph)%state0(:,en) - subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) subFrac = 0.0_pReal - subStep = 1.0_pReal/num%subStepSizeCryst todo = .true. - converged_ = .false. ! pretend failed step of 1/subStepSizeCryst + subStep = 1.0_pReal/num%subStepSizeCryst + converged_ = .false. ! pretend failed step of 1/subStepSizeCryst todo = .true. cutbackLooping: do while (todo) @@ -1038,9 +1033,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en) subState0 = plasticState(ph)%state(:,en) - if (damageState(ph)%sizeState > 0) & - damageState(ph)%subState0(:,en) = damageState(ph)%state(:,en) - endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) @@ -1048,16 +1040,13 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subStep = num%subStepSizeCryst * subStep phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0 phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0 - phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use? + phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use? if (subStep < 1.0_pReal) then ! actual (not initial) cutback phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0 phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0 endif plasticState(ph)%state(:,en) = subState0 - if (damageState(ph)%sizeState > 0) & - damageState(ph)%state(:,en) = damageState(ph)%subState0(:,en) - - todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) + todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif !-------------------------------------------------------------------------------------------------- @@ -1065,13 +1054,12 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) if (todo) then 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,(el-1)*discretization_nIPs + ip) + converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el) endif enddo cutbackLooping -end function crystallite_stress +end function phase_mechanical_constitutive !-------------------------------------------------------------------------------------------------- @@ -1104,12 +1092,13 @@ module subroutine mechanical_restore(ce,includeL) end subroutine mechanical_restore + !-------------------------------------------------------------------------------------------------- !> @brief Calculate tangent (dPdF). !-------------------------------------------------------------------------------------------------- -module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) +module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF) - real(pReal), intent(in) :: dt + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & co, & !< counter in constituent loop ce @@ -1159,11 +1148,11 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal do o=1,3; do p=1,3 lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & - + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * dt + + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * Delta_t lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & + invFi*invFi(p,o) rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & - - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * dt + - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * Delta_t enddo; enddo call math_invert(temp_99,error,math_3333to99(lhs_3333)) if (error) then @@ -1192,7 +1181,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) enddo; enddo - lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * dt & + lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t & + math_mul3333xx3333(dSdFi,dFidS) call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) @@ -1208,7 +1197,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) ! calculate dFpinvdF temp_3333 = math_mul3333xx3333(dLpdS,dSdF) do o=1,3; do p=1,3 - dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * dt + dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * Delta_t enddo; enddo !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 019838689..14d86661b 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -32,7 +32,7 @@ submodule(phase:mechanical) eigen contains -module subroutine eigendeformation_init(phases) +module subroutine eigen_init(phases) class(tNode), pointer :: & phases @@ -68,7 +68,7 @@ module subroutine eigendeformation_init(phases) where(damage_anisobrittle_init()) model_damage = KINEMATICS_cleavage_opening_ID -end subroutine eigendeformation_init +end subroutine eigen_init !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 8d2915ef7..895246dfa 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -216,7 +216,7 @@ module function phase_K_T(co,ce) result(K) end function phase_K_T -module function thermal_stress(Delta_t,ph,en) result(converged_) ! ?? why is this called "stress" when it seems closer to "updateState" ?? +module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) real(pReal), intent(in) :: Delta_t integer, intent(in) :: ph, en @@ -225,7 +225,7 @@ module function thermal_stress(Delta_t,ph,en) result(converged_) ! ?? converged_ = .not. integrateThermalState(Delta_t,ph,en) -end function thermal_stress +end function phase_thermal_constitutive !-------------------------------------------------------------------------------------------------- diff --git a/src/prec.f90 b/src/prec.f90 index 7613cfe46..620647b82 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -38,16 +38,14 @@ module prec sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates offsetDeltaState = 0, & !< index offset of delta state sizeDeltaState = 0 !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments - ! http://stackoverflow.com/questions/3948210 - real(pReal), pointer, dimension(:), contiguous :: & + real(pReal), allocatable, dimension(:) :: & atol - real(pReal), pointer, dimension(:,:), contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous + ! http://stackoverflow.com/questions/3948210 + real(pReal), pointer, dimension(:,:), contiguous :: & !< is basically an allocatable+target, but in a type needs to be pointer state0, & state, & !< state dotState, & !< rate of state change deltaState !< increment of state change - real(pReal), allocatable, dimension(:,:) :: & - subState0 end type type, extends(tState) :: tPlasticState @@ -61,12 +59,9 @@ module prec real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal) !< minimum positive number such that 1.0 + EPSILON /= 1.0. real(pReal), private, parameter :: PREAL_MIN = tiny(0.0_pReal) !< smallest normalized floating point number - integer, dimension(0), parameter :: & - emptyIntArray = [integer::] - real(pReal), dimension(0), parameter :: & - emptyRealArray = [real(pReal)::] - character(len=pStringLen), dimension(0), parameter :: & - emptyStringArray = [character(len=pStringLen)::] + integer, dimension(0), parameter :: emptyIntArray = [integer::] + real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] + character(len=pStringLen), dimension(0), parameter :: emptyStringArray = [character(len=pStringLen)::] private :: & selfTest