From 3f0eafd640ddd60c5009137522ee2a89ef14b164 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 17:53:11 +0200 Subject: [PATCH 01/16] first step towards separating of mechanics, thermal, and damage --- src/CPFEM.f90 | 5 ++-- src/grid/spectral_utilities.f90 | 1 + src/homogenization.f90 | 40 ++++++++++++++++++++++-------- src/mesh/FEM_utilities.f90 | 44 ++++++++++++++++----------------- 4 files changed, 56 insertions(+), 34 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 213b56a54..3a1340764 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -189,9 +189,10 @@ 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 + if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP]) + call materialpoint_stressAndItsTangent2(dt,[ip,ip],[elCP,elCP]) + terminalIllness: if (terminallyIll) then diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 061cc9d74..c143f480f 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -816,6 +816,7 @@ 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 materialpoint_stressAndItsTangent2(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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index f34834272..07218e50a 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -179,6 +179,7 @@ module homogenization public :: & homogenization_init, & materialpoint_stressAndItsTangent, & + materialpoint_stressAndItsTangent2, & homogenization_mu_T, & homogenization_K_T, & homogenization_f_T, & @@ -227,7 +228,7 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- -!> @brief parallelized calculation of stress and corresponding tangent at material points +!> @brief !-------------------------------------------------------------------------------------------------- subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem) @@ -243,8 +244,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE 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) @@ -285,10 +286,30 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE endif enddo enddo - !$OMP END DO + !$OMP END PARALLEL DO + +end subroutine materialpoint_stressAndItsTangent + + +!-------------------------------------------------------------------------------------------------- +!> @brief +!-------------------------------------------------------------------------------------------------- +subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_execElem) + + real(pReal), intent(in) :: dt !< 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 + logical :: & + converged + logical, dimension(2) :: & + doneAndHappy if (.not. terminallyIll) then - !$OMP DO PRIVATE(ho,ph,ce) + !$OMP PARALLEL DO PRIVATE(ho,ph,ce) do el = FEsolving_execElem(1),FEsolving_execElem(2) if (terminallyIll) continue do ip = FEsolving_execIP(1),FEsolving_execIP(2) @@ -305,9 +326,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE enddo enddo enddo - !$OMP END DO + !$OMP END PARALLEL DO - !$OMP DO PRIVATE(ho,ce) + !$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 @@ -318,13 +339,12 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call mechanical_homogenize(dt,ce) enddo IpLooping3 enddo elementLooping3 - !$OMP END DO + !$OMP END PARALLEL DO else print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' endif - !$OMP END PARALLEL -end subroutine materialpoint_stressAndItsTangent +end subroutine materialpoint_stressAndItsTangent2 !-------------------------------------------------------------------------------------------------- diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 6765d3d0d..b388bc95b 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,20 @@ 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 materialpoint_stressAndItsTangent2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) 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 +198,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 From 5f78f1753c7d0e7d138d895c704c3a32577c0df8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 18:03:38 +0200 Subject: [PATCH 02/16] split up thermal only for grid at the moment --- src/grid/spectral_utilities.f90 | 1 + src/homogenization.f90 | 26 +++++++++++++++++++++++++- 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index c143f480f..d164c33a2 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -816,6 +816,7 @@ 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 materialpoint_stressAndItsTangent3(timeinc,[1,1],[1,product(grid(1:2))*grid3]) call materialpoint_stressAndItsTangent2(timeinc,[1,1],[1,product(grid(1:2))*grid3]) P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 07218e50a..9e80c6dda 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -180,6 +180,7 @@ module homogenization homogenization_init, & materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent2, & + materialpoint_stressAndItsTangent3, & homogenization_mu_T, & homogenization_K_T, & homogenization_f_T, & @@ -294,7 +295,7 @@ end subroutine materialpoint_stressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_execElem) +subroutine materialpoint_stressAndItsTangent3(dt,FEsolving_execIP,FEsolving_execElem) real(pReal), intent(in) :: dt !< time increment integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP @@ -327,7 +328,30 @@ subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_exec enddo enddo !$OMP END PARALLEL DO + else + print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' + endif +end subroutine materialpoint_stressAndItsTangent3 +!-------------------------------------------------------------------------------------------------- +!> @brief +!-------------------------------------------------------------------------------------------------- +subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_execElem) + + real(pReal), intent(in) :: dt !< 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 + logical :: & + converged + logical, dimension(2) :: & + doneAndHappy + + + if (.not. terminallyIll) then !$OMP PARALLEL DO PRIVATE(ho,ce) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) From 2a84aa7ae4a612e81d701025e0f55946c24a032c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 20:32:21 +0200 Subject: [PATCH 03/16] obvious, no need for comment --- src/grid/grid_mech_FEM.f90 | 6 ++---- src/grid/grid_mech_spectral_basic.f90 | 6 ++---- src/grid/grid_mech_spectral_polarisation.f90 | 16 +++++++--------- 3 files changed, 11 insertions(+), 17 deletions(-) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 77678137d..d2fd89b2b 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -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..3ec14d00e 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -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..7c9c85199 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -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) From ed6b1be35210087d7da5eaac4b81da22aca69477 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 20:43:08 +0200 Subject: [PATCH 04/16] solver handles terminally ill --- src/CPFEM.f90 | 3 +- src/grid/spectral_utilities.f90 | 8 ++-- src/homogenization.f90 | 67 +++++++++++++++------------------ src/mesh/FEM_utilities.f90 | 5 ++- 4 files changed, 41 insertions(+), 42 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 3a1340764..06d1ca6a8 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -191,7 +191,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS 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]) - call materialpoint_stressAndItsTangent2(dt,[ip,ip],[elCP,elCP]) + if (.not. terminallyIll) & + call materialpoint_stressAndItsTangent2(dt,[ip,ip],[elCP,elCP]) terminalIllness: if (terminallyIll) then diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index d164c33a2..607928c1b 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -816,11 +816,13 @@ 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 materialpoint_stressAndItsTangent3(timeinc,[1,1],[1,product(grid(1:2))*grid3]) - call materialpoint_stressAndItsTangent2(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + if (.not. terminallyIll) & + call materialpoint_stressAndItsTangent3(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + if (.not. terminallyIll) & + call materialpoint_stressAndItsTangent2(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 9e80c6dda..f4f701d7c 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -309,30 +309,28 @@ subroutine materialpoint_stressAndItsTangent3(dt,FEsolving_execIP,FEsolving_exec logical, dimension(2) :: & doneAndHappy - if (.not. terminallyIll) then - !$OMP PARALLEL 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 + !$OMP PARALLEL 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 enddo - !$OMP END PARALLEL DO - else - print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' - endif + enddo + !$OMP END PARALLEL DO + end subroutine materialpoint_stressAndItsTangent3 + !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- @@ -351,22 +349,19 @@ subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_exec doneAndHappy - if (.not. terminallyIll) then - !$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(dt,ce) - enddo IpLooping3 - enddo elementLooping3 - !$OMP END PARALLEL DO - else - print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' - endif + !$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(dt,ce) + enddo IpLooping3 + enddo elementLooping3 + !$OMP END PARALLEL DO + end subroutine materialpoint_stressAndItsTangent2 diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index b388bc95b..402c6ddf7 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -163,8 +163,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) print'(/,a)', ' ... evaluating constitutive response ......................................' call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field - call materialpoint_stressAndItsTangent2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) - cutBack = .false. ! reset cutBack status + if (.not. terminallyIll) & + call materialpoint_stressAndItsTangent2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) + cutBack = .false. 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) From 2b24224c7ee29b04159040e96644a33c9c7a1e28 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 20:50:08 +0200 Subject: [PATCH 05/16] more meaningful name --- src/homogenization.f90 | 4 ++-- src/phase.f90 | 12 ++++++------ src/phase_mechanical.f90 | 4 ++-- src/phase_thermal.f90 | 4 ++-- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index f4f701d7c..57b3dae4b 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -270,7 +270,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) converged = .true. do co = 1, homogenization_Nconstituents(ho) - converged = converged .and. crystallite_stress(dt,co,ip,el) + converged = converged .and. phase_mechanical_constitutive(dt,co,ip,el) enddo if (converged) then @@ -318,7 +318,7 @@ subroutine materialpoint_stressAndItsTangent3(dt,FEsolving_execIP,FEsolving_exec 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. phase_thermal_constitutive(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 diff --git a/src/phase.f90 b/src/phase.f90 index e396e5661..a96145995 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -209,13 +209,13 @@ 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 @@ -225,11 +225,11 @@ module phase logical :: broken end function integrateDamageState - module function crystallite_stress(dt,co,ip,el) result(converged_) + module function phase_mechanical_constitutive(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt integer, intent(in) :: co, ip, el logical :: converged_ - end function crystallite_stress + end function phase_mechanical_constitutive !ToDo: Try to merge the all stiffness functions module function phase_homogenizedC(ph,en) result(C) @@ -315,8 +315,8 @@ module phase plastic_nonlocal_updateCompatibility, & converged, & crystallite_init, & - crystallite_stress, & - thermal_stress, & + phase_mechanical_constitutive, & + phase_thermal_constitutive, & phase_mechanical_dPdF, & crystallite_orientations, & crystallite_push33ToRef, & diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 7562c055f..c176de069 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -977,7 +977,7 @@ end subroutine mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function crystallite_stress(dt,co,ip,el) result(converged_) +module function phase_mechanical_constitutive(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt integer, intent(in) :: & @@ -1071,7 +1071,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) enddo cutbackLooping -end function crystallite_stress +end function phase_mechanical_constitutive !-------------------------------------------------------------------------------------------------- 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 !-------------------------------------------------------------------------------------------------- From 594ad2c310f15fd0676cd5c15d72ec786da07b78 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 22:22:41 +0200 Subject: [PATCH 06/16] semi-separate handling of damage --- src/phase_mechanical.f90 | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index c176de069..2a693c690 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1009,17 +1009,13 @@ module function phase_mechanical_constitutive(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 +1034,6 @@ module function phase_mechanical_constitutive(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 +1041,13 @@ module function phase_mechanical_constitutive(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 !-------------------------------------------------------------------------------------------------- @@ -1066,11 +1056,13 @@ module function phase_mechanical_constitutive(dt,co,ip,el) result(converged_) subF = subF0 & + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) - converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,(el-1)*discretization_nIPs + ip) endif enddo cutbackLooping + if (damageState(ph)%sizeState > 0) damageState(ph)%subState0(:,en) = damageState(ph)%state0(:,en) + converged_ = converged_ .and. .not. integrateDamageState(dt,co,(el-1)*discretization_nIPs + ip) + end function phase_mechanical_constitutive From 777620b800d9de33d4651555b919182496dd43d6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 22:41:38 +0200 Subject: [PATCH 07/16] polishing --- src/phase_damage_anisobrittle.f90 | 14 ++++++-------- src/phase_damage_isobrittle.f90 | 18 ++++++++---------- 2 files changed, 14 insertions(+), 18 deletions(-) diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 9984e5514..e71c6d3cc 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -58,10 +58,10 @@ module function anisobrittle_init() result(mySources) do p = 1, phases%length if(mySources(p)) then - phase => phases%get(p) - sources => phase%get('damage') + phase => phases%get(p) + sources => phase%get('damage') - associate(prm => param(p)) + associate(prm => param(p)) src => sources%get(1) N_cl = src%get_as1dInt('N_cl',defaultVal=emptyIntArray) @@ -97,12 +97,10 @@ module function anisobrittle_init() result(mySources) 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' - 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 diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index aed6a91e9..aad03daa9 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -47,10 +47,10 @@ module function isobrittle_init() result(mySources) 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)) src => sources%get(1) prm%W_crit = src%get_asFloat('W_crit') @@ -69,12 +69,11 @@ 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 + end associate -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)') - endif + + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)') + endif enddo @@ -102,8 +101,7 @@ module subroutine isobrittle_deltaState(C, Fe, ph,me) strain = 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 + strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), & damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me), & From fc735d639190ca09a6cd2413484eff77c54ec5c7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 23:19:48 +0200 Subject: [PATCH 08/16] substate0 only needed for staggered state-stress integration --- src/phase.f90 | 4 ---- src/phase_damage.f90 | 10 +++++----- src/phase_damage_isobrittle.f90 | 5 ++--- src/phase_mechanical.f90 | 4 +--- src/prec.f90 | 17 ++++++----------- 5 files changed, 14 insertions(+), 26 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index a96145995..311e32a1e 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -545,10 +545,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..4cc361637 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -212,8 +212,8 @@ module function integrateDamageState(dt,co,ce) result(broken) 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 + damageState(ph)%state(1:size_so,me) = damageState(ph)%state0 (1:size_so,me) & + + damageState(ph)%dotState(1:size_so,me) * dt source_dotState(1:size_so,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState @@ -228,9 +228,9 @@ module function integrateDamageState(dt,co,ce) result(broken) 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 & + 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 + r(1:size_so) = damageState(ph)%state (1:size_so,me) & + - damageState(ph)%State0 (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) converged_ = converged_ .and. converged(r(1:size_so), & damageState(ph)%state(1:size_so,me), & diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index aad03daa9..6480c3c99 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -103,9 +103,8 @@ module subroutine isobrittle_deltaState(C, Fe, ph,me) associate(prm => param(ph)) strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit - 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)) + damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), 0.0_pReal, & + strainenergy > damageState(ph)%state(1,me)) end associate end subroutine isobrittle_deltaState diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 2a693c690..65dd2ceac 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -197,8 +197,7 @@ module subroutine mechanical_init(materials,phases) Nmembers class(tNode), pointer :: & num_crystallite, & - material, & - constituents, & + phase, & mech @@ -1060,7 +1059,6 @@ module function phase_mechanical_constitutive(dt,co,ip,el) result(converged_) enddo cutbackLooping - if (damageState(ph)%sizeState > 0) damageState(ph)%subState0(:,en) = damageState(ph)%state0(:,en) converged_ = converged_ .and. .not. integrateDamageState(dt,co,(el-1)*discretization_nIPs + ip) end function phase_mechanical_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 From e6d25294d3b1a12e6deef39e8fb7098de0576afd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 23:32:08 +0200 Subject: [PATCH 09/16] separating functionality --- src/phase.f90 | 14 +++++--------- src/phase_damage.f90 | 21 ++++++++++++++++++++- src/phase_damage_isobrittle.f90 | 5 ++--- 3 files changed, 27 insertions(+), 13 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 311e32a1e..9324225a2 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -108,6 +108,10 @@ 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 @@ -435,17 +439,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 diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 4cc361637..700bb2142 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -144,6 +144,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,7 +193,6 @@ 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 diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 6480c3c99..a73fe4fb3 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -121,12 +121,11 @@ module subroutine isobrittle_results(phase,group) integer :: o - associate(prm => param(phase), & - stt => damageState(phase)%state) + associate(prm => param(phase), stt => damageState(phase)%state) 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 From f9edeb40a5a3e3fb2a15584365ab16c90022c91f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 17 Jul 2021 11:50:21 +0200 Subject: [PATCH 10/16] descriptive names --- src/CPFEM.f90 | 4 +- src/grid/grid_mech_FEM.f90 | 2 +- src/grid/grid_mech_spectral_basic.f90 | 2 +- src/grid/grid_mech_spectral_polarisation.f90 | 2 +- src/grid/spectral_utilities.f90 | 6 +- src/homogenization.f90 | 74 ++++++++------------ src/homogenization_mechanical.f90 | 8 +-- src/mesh/FEM_utilities.f90 | 4 +- src/phase.f90 | 14 ++-- src/phase_damage.f90 | 14 ++-- src/phase_damage_anisobrittle.f90 | 18 ++--- src/phase_mechanical.f90 | 27 +++---- src/phase_mechanical_eigen.f90 | 4 +- 13 files changed, 83 insertions(+), 96 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 06d1ca6a8..aa532859a 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -190,9 +190,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS 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]) + call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP]) if (.not. terminallyIll) & - call materialpoint_stressAndItsTangent2(dt,[ip,ip],[elCP,elCP]) + 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 d2fd89b2b..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 diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 3ec14d00e..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 diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 7c9c85199..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 diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 607928c1b..bd37e2654 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -815,11 +815,11 @@ 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 materialpoint_stressAndItsTangent3(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + call homogenization_thermal_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) if (.not. terminallyIll) & - call materialpoint_stressAndItsTangent2(timeinc,[1,1],[1,product(grid(1:2))*grid3]) + 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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 57b3dae4b..03a7c620d 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,9 +178,9 @@ module homogenization public :: & homogenization_init, & - materialpoint_stressAndItsTangent, & - materialpoint_stressAndItsTangent2, & - materialpoint_stressAndItsTangent3, & + homogenization_mechanical_response, & + homogenization_mechanical_response2, & + homogenization_thermal_response, & homogenization_mu_T, & homogenization_K_T, & homogenization_f_T, & @@ -231,15 +231,15 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- !> @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) :: & @@ -268,13 +268,10 @@ 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. phase_mechanical_constitutive(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.] @@ -282,34 +279,30 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE enddo convergenceLooping 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 PARALLEL DO -end subroutine materialpoint_stressAndItsTangent +end subroutine homogenization_mechanical_response !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine materialpoint_stressAndItsTangent3(dt,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_thermal_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 - logical :: & - converged - logical, dimension(2) :: & - doneAndHappy + ip, & !< integration point number + el, & !< element number + co, ce, ho - !$OMP PARALLEL DO PRIVATE(ho,ph,ce) + + !$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) @@ -317,36 +310,29 @@ subroutine materialpoint_stressAndItsTangent3(dt,FEsolving_execIP,FEsolving_exec ho = material_homogenizationID(ce) call thermal_partition(ce) do co = 1, homogenization_Nconstituents(ho) - ph = material_phaseID(co,ce) - if (.not. phase_thermal_constitutive(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 + 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 enddo !$OMP END PARALLEL DO -end subroutine materialpoint_stressAndItsTangent3 +end subroutine homogenization_thermal_response !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_execElem) +subroutine homogenization_mechanical_response2(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 - logical :: & - converged - logical, dimension(2) :: & - doneAndHappy + ip, & !< integration point number + el, & !< element number + co, ce, ho !$OMP PARALLEL DO PRIVATE(ho,ce) @@ -357,13 +343,13 @@ subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_exec do co = 1, homogenization_Nconstituents(ho) call crystallite_orientations(co,ip,el) enddo - call mechanical_homogenize(dt,ce) + call mechanical_homogenize(Delta_t,ce) enddo IpLooping3 enddo elementLooping3 !$OMP END PARALLEL DO -end subroutine materialpoint_stressAndItsTangent2 +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 402c6ddf7..5bb6d5f44 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -162,9 +162,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) 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 materialpoint_stressAndItsTangent2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) + call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) cutBack = .false. P_av = sum(homogenization_P,dim=3) * wgt diff --git a/src/phase.f90 b/src/phase.f90 index 9324225a2..b551c4170 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -113,8 +113,8 @@ module phase 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 @@ -221,21 +221,21 @@ module phase end function phase_thermal_constitutive - module function integrateDamageState(dt,co,ce) result(broken) - real(pReal), intent(in) :: dt + module function integrateDamageState(Delta_t,co,ce) result(broken) + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & ce, & co logical :: broken end function integrateDamageState - module function phase_mechanical_constitutive(dt,co,ip,el) result(converged_) - real(pReal), intent(in) :: dt + 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: Try to merge the all stiffness functions + !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 diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 700bb2142..89d3b1653 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -197,9 +197,9 @@ 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) +module function integrateDamageState(Delta_t,co,ce) result(broken) - real(pReal), intent(in) :: dt + real(pReal), intent(in) :: Delta_t integer, intent(in) :: & ce, & co @@ -230,10 +230,10 @@ module function integrateDamageState(dt,co,ce) result(broken) broken = phase_damage_collectDotState(ph,me) if(broken) return - size_so = damageState(ph)%sizeDotState - damageState(ph)%state(1:size_so,me) = damageState(ph)%state0 (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,me) = damageState(ph)%state0 (1:size_so,me) & + + damageState(ph)%dotState(1:size_so,me) * Delta_t + source_dotState(1:size_so,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState @@ -249,7 +249,7 @@ module function integrateDamageState(dt,co,ce) result(broken) + source_dotState(1:size_so,1)* (1.0_pReal - zeta) r(1:size_so) = damageState(ph)%state (1:size_so,me) & - damageState(ph)%State0 (1:size_so,me) & - - damageState(ph)%dotState(1:size_so,me) * dt + - damageState(ph)%dotState(1:size_so,me) * Delta_t damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so) converged_ = converged_ .and. converged(r(1:size_so), & damageState(ph)%state(1:size_so,me), & diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index e71c6d3cc..cc3628323 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) + 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,10 +92,10 @@ 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 diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 65dd2ceac..e03998989 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 @@ -302,7 +302,7 @@ module subroutine mechanical_init(materials,phases) end select - call eigendeformation_init(phases) + call eigen_init(phases) end subroutine mechanical_init @@ -976,9 +976,9 @@ end subroutine mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function phase_mechanical_constitutive(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, & @@ -1054,12 +1054,12 @@ module function phase_mechanical_constitutive(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_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el) endif enddo cutbackLooping - converged_ = converged_ .and. .not. integrateDamageState(dt,co,(el-1)*discretization_nIPs + ip) + converged_ = converged_ .and. .not. integrateDamageState(Delta_t,co,(el-1)*discretization_nIPs + ip) end function phase_mechanical_constitutive @@ -1094,12 +1094,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 @@ -1149,11 +1150,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 @@ -1182,7 +1183,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)) @@ -1198,7 +1199,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 !-------------------------------------------------------------------------------------------------- From c109d5a37b349e937c92c09901264471ef4f28db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 17 Jul 2021 14:06:48 +0200 Subject: [PATCH 11/16] better have different physics separated --- src/homogenization.f90 | 4 ++-- src/phase.f90 | 12 +++++------- src/phase_damage.f90 | 23 ++++++++++++++++++++++- src/phase_mechanical.f90 | 2 -- 4 files changed, 29 insertions(+), 12 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 03a7c620d..12d192025 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -233,7 +233,7 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem) - real(pReal), intent(in) :: Delta_t !< time increment + real(pReal), intent(in) :: Delta_t !< time increment integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP integer :: & NiterationMPstate, & @@ -269,7 +269,7 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) - + converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) if (converged) then doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) converged = all(doneAndHappy) diff --git a/src/phase.f90 b/src/phase.f90 index b551c4170..49f1adc16 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -221,13 +221,11 @@ module phase end function phase_thermal_constitutive - module function integrateDamageState(Delta_t,co,ce) result(broken) + module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) real(pReal), intent(in) :: Delta_t - integer, intent(in) :: & - ce, & - co - logical :: broken - end function integrateDamageState + integer, intent(in) :: co, ip, el + logical :: converged_ + end function phase_damage_constitutive module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_) real(pReal), intent(in) :: Delta_t @@ -321,12 +319,12 @@ module phase crystallite_init, & 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, & diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 89d3b1653..4164854fe 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -125,6 +125,27 @@ 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 + + + converged_ = .not. integrateDamageState(Delta_t,co,(el-1)*discretization_nIPs + ip) + +end function phase_damage_constitutive + + !-------------------------------------------------------------------------------------------------- !> @brief returns the degraded/modified elasticity matrix !-------------------------------------------------------------------------------------------------- @@ -197,7 +218,7 @@ 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(Delta_t,co,ce) result(broken) +function integrateDamageState(Delta_t,co,ce) result(broken) real(pReal), intent(in) :: Delta_t integer, intent(in) :: & diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index e03998989..8b04a0009 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1059,8 +1059,6 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged enddo cutbackLooping - converged_ = converged_ .and. .not. integrateDamageState(Delta_t,co,(el-1)*discretization_nIPs + ip) - end function phase_mechanical_constitutive From d068f45aa097d6e0be83a25641d57686a356613b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 17 Jul 2021 15:25:00 +0200 Subject: [PATCH 12/16] avoid superflouos damage calculations --- src/homogenization.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 12d192025..fc5dd7bda 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -269,15 +269,16 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) - converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))]) if (converged) then 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*, ' Cell ', ce, ' terminally ill' terminallyIll = .true. From 50fccb0a2ea7e0fe893650ac9d92ae7c84797ea8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 17 Jul 2021 22:53:58 +0200 Subject: [PATCH 13/16] not needed --- examples/config/phase/damage/isobrittle_generic.yaml | 2 -- 1 file changed, 2 deletions(-) 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] From 57ad308a7f46cfd91360bc3f67bdff706064822b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 18 Jul 2021 09:10:33 +0200 Subject: [PATCH 14/16] readable --- src/phase_damage_isobrittle.f90 | 35 +++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index a73fe4fb3..1c3196a15 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') - 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,6 +79,9 @@ 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' + stt%r_W => damageState(ph)%state(1,:) + dlt%r_W => damageState(ph)%deltaState(1,:) + end associate @@ -93,18 +106,18 @@ module subroutine isobrittle_deltaState(C, Fe, ph,me) 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*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(me) = merge(r_W - stt%r_W(me), 0.0_pReal, r_W > stt%r_W(me)) - damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), 0.0_pReal, & - strainenergy > damageState(ph)%state(1,me)) end associate end subroutine isobrittle_deltaState @@ -121,13 +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³') ! Wrong, this is dimensionless end select enddo outputsLoop + end associate end subroutine isobrittle_results From 6ad6158bfb590830321b4f6e8afc583948741d0d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 18 Jul 2021 09:44:52 +0200 Subject: [PATCH 15/16] (en)try is the name used in the DADF5 file --- src/phase.f90 | 8 ++-- src/phase_damage.f90 | 64 +++++++++++++++---------------- src/phase_damage_anisobrittle.f90 | 16 ++++---- src/phase_damage_isobrittle.f90 | 6 +-- 4 files changed, 47 insertions(+), 47 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 49f1adc16..d90f41ccc 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -168,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 @@ -273,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) :: & diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 4164854fe..12a5ddd41 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 @@ -229,7 +229,7 @@ function integrateDamageState(Delta_t,co,ce) result(broken) integer :: & NiterationState, & !< number of iterations in state loop ph, & - me, & + en, & size_so real(pReal) :: & zeta @@ -240,7 +240,7 @@ function integrateDamageState(Delta_t,co,ce) result(broken) converged_ ph = material_phaseID(co,ce) - me = material_phaseEntry(co,ce) + en = material_phaseEntry(co,ce) if (damageState(ph)%sizeState == 0) then broken = .false. @@ -248,37 +248,37 @@ function integrateDamageState(Delta_t,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)%state0 (1:size_so,me) & - + damageState(ph)%dotState(1:size_so,me) * Delta_t + 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)%State0 (1:size_so,me) & - - damageState(ph)%dotState(1:size_so,me) * Delta_t - 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 @@ -340,11 +340,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 @@ -355,11 +355,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 @@ -398,11 +398,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 :: & @@ -419,13 +419,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 @@ -476,13 +476,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 cc3628323..213d04f28 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -110,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 @@ -124,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 + & @@ -169,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) :: & @@ -191,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 1c3196a15..066cab47e 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -97,9 +97,9 @@ 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) :: & @@ -116,7 +116,7 @@ module subroutine isobrittle_deltaState(C, Fe, ph,me) 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(me) = merge(r_W - stt%r_W(me), 0.0_pReal, r_W > stt%r_W(me)) + dlt%r_W(en) = merge(r_W - stt%r_W(en), 0.0_pReal, r_W > stt%r_W(en)) end associate From 53b7aab29dc3ba165a65e4038bd71f6eb7818faa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 18 Jul 2021 09:48:49 +0200 Subject: [PATCH 16/16] use ph,en access pattern --- src/phase_damage.f90 | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 12a5ddd41..172f6f22a 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -140,8 +140,10 @@ module function phase_damage_constitutive(Delta_t,co,ip,el) result(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,co,(el-1)*discretization_nIPs + ip) + converged_ = .not. integrateDamageState(Delta_t,ph,en) end function phase_damage_constitutive @@ -218,18 +220,16 @@ end function phase_f_phi !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -function integrateDamageState(Delta_t,co,ce) result(broken) +function integrateDamageState(Delta_t,ph,en) result(broken) real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - ce, & - co + ph, & + en logical :: broken integer :: & NiterationState, & !< number of iterations in state loop - ph, & - en, & size_so real(pReal) :: & zeta @@ -239,8 +239,6 @@ function integrateDamageState(Delta_t,co,ce) result(broken) logical :: & converged_ - ph = material_phaseID(co,ce) - en = material_phaseEntry(co,ce) if (damageState(ph)%sizeState == 0) then broken = .false.