From 3f0eafd640ddd60c5009137522ee2a89ef14b164 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 16 Jul 2021 17:53:11 +0200 Subject: [PATCH] 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