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

store field variables as 1D array 

first step of simplifying layout:
1) Solver translates from ip,el tuple (FEM) or cells(1),cells(2),cells(3) triple to list.
2) DAMASK iterates over all points
3) homogenization knows mapping (point,constituent) -> (instance,member)
This commit is contained in:
Martin Diehl 2020-12-16 12:48:45 +01:00
parent 5d9c931008
commit 3884549e19
10 changed files with 86 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/CPFEM.f90
View File

@ -153,7 +153,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
H H
integer(pInt) elCP, & ! crystal plasticity element number integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog, mySource i, j, k, l, m, n, ph, homog, mySource,ma
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll
@ -161,6 +161,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
elCP = mesh_FEM2DAMASK_elem(elFE) elCP = mesh_FEM2DAMASK_elem(elFE)
ma = (elCP-1) * discretization_nIPs + ip
if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then
print'(/,a)', '#############################################' print'(/,a)', '#############################################'
print'(a1,a22,1x,i8,a13)', '#','element', elCP, '#' print'(a1,a22,1x,i8,a13)', '#','element', elCP, '#'
@ -184,8 +186,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
temperature_inp temperature_inp
end select chosenThermal1 end select chosenThermal1
homogenization_F0(1:3,1:3,ip,elCP) = ffn homogenization_F0(1:3,1:3,ma) = ffn
homogenization_F(1:3,1:3,ip,elCP) = ffn1 homogenization_F(1:3,1:3,ma) = ffn1
if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then
@ -212,17 +214,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else terminalIllness else terminalIllness
! translate from P to sigma ! translate from P to sigma
Kirchhoff = matmul(homogenization_P(1:3,1:3,ip,elCP), transpose(homogenization_F(1:3,1:3,ip,elCP))) Kirchhoff = matmul(homogenization_P(1:3,1:3,ma), transpose(homogenization_F(1:3,1:3,ma)))
J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ip,elCP)) J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ma))
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
! translate from dP/dF to dCS/dE ! translate from dP/dF to dCS/dE
H = 0.0_pReal H = 0.0_pReal
do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
H(i,j,k,l) = H(i,j,k,l) & H(i,j,k,l) = H(i,j,k,l) &
+ homogenization_F(j,m,ip,elCP) * homogenization_F(l,n,ip,elCP) & + homogenization_F(j,m,ma) * homogenization_F(l,n,ma) &
* homogenization_dPdF(i,m,k,n,ip,elCP) & * homogenization_dPdF(i,m,k,n,ma) &
- math_delta(j,l) * homogenization_F(i,m,ip,elCP) * homogenization_P(k,m,ip,elCP) & - math_delta(j,l) * homogenization_F(i,m,ma) * homogenization_P(k,m,ma) &
+ 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) & + 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
+ Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k)) + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo

30
src/grid/grid_mech_FEM.f90
View File

@ -238,7 +238,7 @@ subroutine grid_mech_FEM_init
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restartRead endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,1,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 materialpoint_stressAndItsTangent
call utilities_updateCoords(F) 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 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 F, & ! target F
@ -359,7 +359,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
F_lastInc = F F_lastInc = F
homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F, [3,3,product(grid(1:2))*grid3])
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -557,9 +557,9 @@ subroutine formResidual(da_local,x_local, &
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
ele = ele + 1 ele = ele + 1
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + &
matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,1,ele) + & matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
ctr = 0 ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
ctr = ctr + 1 ctr = ctr + 1
@ -636,18 +636,18 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
row = col row = col
ele = ele + 1 ele = ele + 1
K_ele = 0.0 K_ele = 0.0
K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + & K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + & K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + & K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
K_ele = K_ele + & K_ele = K_ele + &
matmul(transpose(BMatFull), & matmul(transpose(BMatFull), &
matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,1,ele), & matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,ele), &
shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ
call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)

4
src/grid/grid_mech_spectral_basic.f90
View File

@ -199,7 +199,7 @@ subroutine grid_mech_spectral_basic_init
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restartRead endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,1,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 materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc))) 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 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 reshape(F,shape(F_lastInc)), & ! target F
@ -319,7 +319,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_
rotation_BC%rotate(F_aimDot,active=.true.)) rotation_BC%rotate(F_aimDot,active=.true.))
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

6
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -225,7 +225,7 @@ subroutine grid_mech_spectral_polarisation_init
F_tau_lastInc = 2.0_pReal*F_lastInc F_tau_lastInc = 2.0_pReal*F_lastInc
endif restartRead endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,1,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 materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc))) 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 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 reshape(F,shape(F_lastInc)), & ! target F
@ -359,7 +359,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -604,7 +604,7 @@ subroutine formResidual(in, FandF_tau, &
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
e = e + 1 e = e + 1
residual_F(1:3,1:3,i,j,k) = & residual_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
+ residual_F_tau(1:3,1:3,i,j,k) + residual_F_tau(1:3,1:3,i,j,k)

16
src/grid/spectral_utilities.f90
View File

@ -810,7 +810,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
print'(/,a)', ' ... evaluating constitutive response ......................................' print'(/,a)', ' ... evaluating constitutive response ......................................'
flush(IO_STDOUT) flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
call materialpoint_stressAndItsTangent(timeinc) ! calculate P field call materialpoint_stressAndItsTangent(timeinc) ! calculate P field
@ -829,13 +829,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
dPdF_min = huge(1.0_pReal) dPdF_min = huge(1.0_pReal)
dPdF_norm_min = huge(1.0_pReal) dPdF_norm_min = huge(1.0_pReal)
do i = 1, product(grid(1:2))*grid3 do i = 1, product(grid(1:2))*grid3
if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) then
dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i) dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)
endif endif
if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) then
dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i) dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)
endif endif
end do end do
@ -853,7 +853,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min) C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min)
C_volAvg = sum(sum(homogenization_dPdF,dim=6),dim=5) C_volAvg = sum(homogenization_dPdF,dim=5)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (ierr /= 0) error stop 'MPI error' if (ierr /= 0) error stop 'MPI error'
C_volAvg = C_volAvg * wgt C_volAvg = C_volAvg * wgt

23
src/homogenization.f90
View File

@ -30,12 +30,12 @@ module homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point ! General variables for the homogenization at a material point
real(pReal), dimension(:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:), allocatable, public :: &
homogenization_F0, & !< def grad of IP at start of FE increment homogenization_F0, & !< def grad of IP at start of FE increment
homogenization_F !< def grad of IP to be reached at end of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment
real(pReal), dimension(:,:,:,:), allocatable, public :: & !, protected :: & ! Issue with ifort real(pReal), dimension(:,:,:), allocatable, public :: & !, protected :: & Issue with ifort
homogenization_P !< first P--K stress of IP homogenization_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: & !, protected :: & real(pReal), dimension(:,:,:,:,:), allocatable, public :: & !, protected :: &
homogenization_dPdF !< tangent of first P--K stress at IP homogenization_dPdF !< tangent of first P--K stress at IP
@ -193,6 +193,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
converged converged
logical, dimension(2,discretization_nIPs,discretization_Nelems) :: & logical, dimension(2,discretization_nIPs,discretization_Nelems) :: &
doneAndHappy doneAndHappy
integer :: m
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -227,7 +228,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),& any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),&
FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog))
!$OMP PARALLEL DO !$OMP PARALLEL DO PRIVATE(m)
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Nconstituents(material_homogenizationAt(e)) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2)
@ -297,13 +298,14 @@ subroutine materialpoint_stressAndItsTangent(dt)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! deformation partitioning ! deformation partitioning
!$OMP PARALLEL DO PRIVATE(myNgrains) !$OMP PARALLEL DO PRIVATE(myNgrains,m)
elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Nconstituents(material_homogenizationAt(e)) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done
call mech_partition(homogenization_F0(1:3,1:3,i,e) & m = (e-1)*discretization_nIPs + i
+ (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e))& call mech_partition(homogenization_F0(1:3,1:3,m) &
+ (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m))&
*(subStep(i,e)+subFrac(i,e)), & *(subStep(i,e)+subFrac(i,e)), &
i,e) i,e)
crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains
@ -321,16 +323,17 @@ subroutine materialpoint_stressAndItsTangent(dt)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state update ! state update
!$OMP PARALLEL DO !$OMP PARALLEL DO PRIVATE(m)
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then
if (.not. converged(i,e)) then if (.not. converged(i,e)) then
doneAndHappy(1:2,i,e) = [.true.,.false.] doneAndHappy(1:2,i,e) = [.true.,.false.]
else else
m = (e-1)*discretization_nIPs + i
doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), & doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), &
homogenization_F0(1:3,1:3,i,e) & homogenization_F0(1:3,1:3,m) &
+ (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e)) & + (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m)) &
*(subStep(i,e)+subFrac(i,e)), & *(subStep(i,e)+subFrac(i,e)), &
i,e) i,e)
converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy

23
src/homogenization_mech.f90
View File

@ -73,10 +73,10 @@ module subroutine mech_init(num_homog)
print'(/,a)', ' <<<+- homogenization_mech init -+>>>' print'(/,a)', ' <<<+- homogenization_mech init -+>>>'
allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
homogenization_F0 = spread(spread(math_I3,3,discretization_nIPs),4,discretization_Nelems) ! initialize to identity homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity
homogenization_F = homogenization_F0 ! initialize to identity homogenization_F = homogenization_F0 ! initialize to identity
allocate(homogenization_P(3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal) allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict) num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init
@ -127,23 +127,24 @@ module subroutine mech_homogenize(ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
integer :: c integer :: c,m
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
m = (el-1)* discretization_nIPs + ip
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
homogenization_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) homogenization_P(1:3,1:3,m) = crystallite_P(1:3,1:3,1,ip,el)
homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el) homogenization_dPdF(1:3,1:3,1:3,1:3,m) = crystallite_stressTangent(1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do c = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo enddo
call mech_isostrain_averageStressAndItsTangent(& call mech_isostrain_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ip,el), & homogenization_P(1:3,1:3,m), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, & dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
@ -153,8 +154,8 @@ module subroutine mech_homogenize(ip,el)
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo enddo
call mech_RGC_averageStressAndItsTangent(& call mech_RGC_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ip,el), & homogenization_P(1:3,1:3,m), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),& homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, & dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))

1
src/material.f90
View File

@ -140,7 +140,6 @@ contains
subroutine material_init(restart) subroutine material_init(restart)
logical, intent(in) :: restart logical, intent(in) :: restart
integer :: myHomog
print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT)

2
src/mesh/FEM_utilities.f90
View File

@ -164,7 +164,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
cutBack = .false. ! reset cutBack status cutBack = .false. ! reset cutBack status
P_av = sum(sum(homogenization_P,dim=4),dim=3) * wgt ! average of P P_av = sum(homogenization_P,dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine utilities_constitutiveResponse end subroutine utilities_constitutiveResponse

43
src/mesh/mesh_mech_FEM.f90
View File

@ -316,16 +316,16 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
Vec :: x_local, f_local, xx_local Vec :: x_local, f_local, xx_local
PetscSection :: section PetscSection :: section
PetscScalar, dimension(:), pointer :: x_scal, pf_scal PetscScalar, dimension(:), pointer :: x_scal, pf_scal
PetscScalar, target :: f_scal(cellDof) PetscScalar, dimension(cellDof), target :: f_scal
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex) PetscReal :: IcellJMat(dimPlex,dimPlex)
PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx, & qPt, basis, comp, cidx, &
numFields numFields, &
PetscReal :: detFAvg bcSize,m
PetscReal :: BMat(dimPlex*dimPlex,cellDof) PetscReal :: detFAvg, detJ
PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat
PetscInt :: bcSize
IS :: bcPoints IS :: bcPoints
@ -366,6 +366,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
do qPt = 0, nQuadrature-1 do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1
BMat = 0.0 BMat = 0.0
do basis = 0, nBasis-1 do basis = 0, nBasis-1
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
@ -375,15 +376,14 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo enddo
if (num%BBarStabilisation) then if (num%BBarStabilisation) then
detFAvg = math_det33(sum(homogenization_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature))
do qPt = 1, nQuadrature do qPt = 0, nQuadrature-1
homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1) = & m = cell*nQuadrature + qPt+1
homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1)* & homogenization_F(1:dimPlex,1:dimPlex,m) = homogenization_F(1:dimPlex,1:dimPlex,m) &
(detFAvg/math_det33(homogenization_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex)) * (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0/real(dimPlex))
enddo enddo
endif endif
@ -407,6 +407,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
f_scal = 0.0 f_scal = 0.0
do qPt = 0, nQuadrature-1 do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1
BMat = 0.0 BMat = 0.0
do basis = 0, nBasis-1 do basis = 0, nBasis-1
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
@ -418,7 +419,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
enddo enddo
f_scal = f_scal + & f_scal = f_scal + &
matmul(transpose(BMat), & matmul(transpose(BMat), &
reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), & reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), &
shape=[dimPlex*dimPlex]))*qWeights(qPt+1) shape=[dimPlex*dimPlex]))*qWeights(qPt+1)
enddo enddo
f_scal = f_scal*abs(detJ) f_scal = f_scal*abs(detJ)
@ -463,7 +464,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
K_eB K_eB
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx,bcSize qPt, basis, comp, cidx,bcSize, m
IS :: bcPoints IS :: bcPoints
@ -506,6 +507,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
FAvg = 0.0 FAvg = 0.0
BMatAvg = 0.0 BMatAvg = 0.0
do qPt = 0, nQuadrature-1 do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt + 1
BMat = 0.0 BMat = 0.0
do basis = 0, nBasis-1 do basis = 0, nBasis-1
do comp = 0, dimPlex-1 do comp = 0, dimPlex-1
@ -516,7 +518,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
if (num%BBarStabilisation) then if (num%BBarStabilisation) then
@ -524,12 +526,11 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
FInv = math_inv33(F) FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
K_eB = K_eB - & K_eB = K_eB - &
matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), & matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex*dimPlex,1]), &
shape=[dimPlex*dimPlex,1]), &
matmul(reshape(FInv(1:dimPlex,1:dimPlex), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA) shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA)
MatB = MatB + & MatB = MatB &
matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA) + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1,dimPlex*dimPlex]),MatA)
FAvg = FAvg + F FAvg = FAvg + F
BMatAvg = BMatAvg + BMat BMatAvg = BMatAvg + BMat
else else