223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

first step towards separating of mechanics, thermal, and damage

This commit is contained in:
Martin Diehl 2021-07-16 17:53:11 +02:00
parent e2b8145dc6
commit 3f0eafd640
4 changed files with 56 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
src/CPFEM.f90
View File

@ -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) CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
else validCalculation else validCalculation
if (debugCPFEM%extensive) & if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP]) call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP])
call materialpoint_stressAndItsTangent2(dt,[ip,ip],[elCP,elCP])
terminalIllness: if (terminallyIll) then terminalIllness: if (terminallyIll) then

1
src/grid/spectral_utilities.f90
View File

@ -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 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_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 = 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 ! average of P

40
src/homogenization.f90
View File

@ -179,6 +179,7 @@ module homogenization
public :: & public :: &
homogenization_init, & homogenization_init, &
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
materialpoint_stressAndItsTangent2, &
homogenization_mu_T, & homogenization_mu_T, &
homogenization_K_T, & homogenization_K_T, &
homogenization_f_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) subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem)
@ -243,8 +244,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
logical, dimension(2) :: & logical, dimension(2) :: &
doneAndHappy 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 el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2) do ip = FEsolving_execIP(1),FEsolving_execIP(2)
@ -285,10 +286,30 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
endif endif
enddo enddo
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 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) do el = FEsolving_execElem(1),FEsolving_execElem(2)
if (terminallyIll) continue if (terminallyIll) continue
do ip = FEsolving_execIP(1),FEsolving_execIP(2) do ip = FEsolving_execIP(1),FEsolving_execIP(2)
@ -305,9 +326,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
enddo enddo
enddo 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) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip ce = (el-1)*discretization_nIPs + ip
@ -318,13 +339,12 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
call mechanical_homogenize(dt,ce) call mechanical_homogenize(dt,ce)
enddo IpLooping3 enddo IpLooping3
enddo elementLooping3 enddo elementLooping3
!$OMP END DO !$OMP END PARALLEL DO
else else
print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' print'(/,a,/)', ' << HOMOG >> Material Point terminally ill'
endif endif
!$OMP END PARALLEL
end subroutine materialpoint_stressAndItsTangent end subroutine materialpoint_stressAndItsTangent2
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

44
src/mesh/FEM_utilities.f90
View File

@ -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 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, parameter :: maxFields = 6
integer, public :: nActiveFields = 0 integer, public :: nActiveFields = 0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! grid related information information ! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems real(pReal), public :: wgt !< weighting factor 1/Nelems
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! field labels information ! field labels information
character(len=*), parameter, public :: & character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical' FIELD_MECH_label = 'mechanical'
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
FIELD_UNDEFINED_ID, & FIELD_UNDEFINED_ID, &
FIELD_MECH_ID FIELD_MECH_ID
@ -48,7 +48,7 @@ module FEM_utilities
COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID COMPONENT_MECH_Z_ID
end enum end enum
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables controlling debugging ! variables controlling debugging
logical :: & logical :: &
@ -57,23 +57,23 @@ module FEM_utilities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type, public :: tSolutionState !< return type of solution from FEM solver variants type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true. logical :: converged = .true.
logical :: stagConverged = .true. logical :: stagConverged = .true.
integer :: iterationsNeeded = 0 integer :: iterationsNeeded = 0
end type tSolutionState end type tSolutionState
type, public :: tComponentBC type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable, dimension(:) :: Value real(pReal), allocatable, dimension(:) :: Value
logical, allocatable, dimension(:) :: Mask logical, allocatable, dimension(:) :: Mask
end type tComponentBC end type tComponentBC
type, public :: tFieldBC type, public :: tFieldBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer :: nComponents = 0 integer :: nComponents = 0
type(tComponentBC), allocatable :: componentBC(:) type(tComponentBC), allocatable :: componentBC(:)
end type tFieldBC end type tFieldBC
type, public :: tLoadCase type, public :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments integer :: incs = 0, & !< number of increments
@ -83,7 +83,7 @@ module FEM_utilities
integer, allocatable, dimension(:) :: faceID integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable, dimension(:) :: fieldBC type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase end type tLoadCase
public :: & public :: &
FEM_utilities_init, & FEM_utilities_init, &
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
@ -94,14 +94,14 @@ module FEM_utilities
COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID COMPONENT_MECH_Z_ID
contains contains
!ToDo: use functions in variable call !ToDo: use functions in variable call
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, sets debug flags !> @brief allocates all neccessary fields, sets debug flags
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_utilities_init subroutine FEM_utilities_init
character(len=pStringLen) :: petsc_optionsOrder character(len=pStringLen) :: petsc_optionsOrder
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh, & num_mesh, &
@ -113,7 +113,7 @@ subroutine FEM_utilities_init
PetscErrorCode :: ierr PetscErrorCode :: ierr
print'(/,a)', ' <<<+- FEM_utilities init -+>>>' print'(/,a)', ' <<<+- FEM_utilities init -+>>>'
num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2) 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 write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal) wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
@ -152,20 +152,20 @@ end subroutine FEM_utilities_init
!> @brief calculates constitutive response !> @brief calculates constitutive response
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
real(pReal), intent(in) :: timeinc !< loading time real(pReal), intent(in) :: timeinc !< loading time
logical, intent(in) :: forwardData !< age results logical, intent(in) :: forwardData !< age results
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
PetscErrorCode :: ierr PetscErrorCode :: ierr
print'(/,a)', ' ... evaluating constitutive response ......................................' print'(/,a)', ' ... evaluating constitutive response ......................................'
call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field 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 cutBack = .false. ! reset cutBack status
P_av = sum(homogenization_P,dim=3) * wgt 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) 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 do dof = offset+comp+1, offset+numDof, numComp
localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
enddo enddo
enddo enddo
call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr) call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr) call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr)
call VecAssemblyEnd (localVec, ierr); CHKERRQ(ierr) call VecAssemblyEnd (localVec, ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr) if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr)
end subroutine utilities_projectBCValues end subroutine utilities_projectBCValues
end module FEM_utilities end module FEM_utilities