solver handles terminally ill

This commit is contained in:
Martin Diehl 2021-07-16 20:43:08 +02:00
parent 2a84aa7ae4
commit ed6b1be352
4 changed files with 41 additions and 42 deletions

View File

@ -191,7 +191,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else validCalculation 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_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 terminalIllness: if (terminallyIll) then

View File

@ -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 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_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]) 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 = 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) 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))', & 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 ' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal

View File

@ -309,30 +309,28 @@ subroutine materialpoint_stressAndItsTangent3(dt,FEsolving_execIP,FEsolving_exec
logical, dimension(2) :: & logical, dimension(2) :: &
doneAndHappy doneAndHappy
if (.not. terminallyIll) then !$OMP PARALLEL 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) ce = (el-1)*discretization_nIPs + ip
ce = (el-1)*discretization_nIPs + ip ho = material_homogenizationID(ce)
ho = material_homogenizationID(ce) call thermal_partition(ce)
call thermal_partition(ce) do co = 1, homogenization_Nconstituents(ho)
do co = 1, homogenization_Nconstituents(ho) ph = material_phaseID(co,ce)
ph = material_phaseID(co,ce) if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then
if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then if (.not. terminallyIll) & ! so first signals terminally ill...
if (.not. terminallyIll) & ! so first signals terminally ill... print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
print*, ' Integration point ', ip,' at element ', el, ' terminally ill' terminallyIll = .true. ! ...and kills all others
terminallyIll = .true. ! ...and kills all others endif
endif
enddo
enddo enddo
enddo enddo
!$OMP END PARALLEL DO enddo
else !$OMP END PARALLEL DO
print'(/,a,/)', ' << HOMOG >> Material Point terminally ill'
endif
end subroutine materialpoint_stressAndItsTangent3 end subroutine materialpoint_stressAndItsTangent3
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief !> @brief
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -351,22 +349,19 @@ subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_exec
doneAndHappy doneAndHappy
if (.not. terminallyIll) then !$OMP PARALLEL 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 ho = material_homogenizationID(ce)
ho = material_homogenizationID(ce) do co = 1, homogenization_Nconstituents(ho)
do co = 1, homogenization_Nconstituents(ho) call crystallite_orientations(co,ip,el)
call crystallite_orientations(co,ip,el) enddo
enddo call mechanical_homogenize(dt,ce)
call mechanical_homogenize(dt,ce) enddo IpLooping3
enddo IpLooping3 enddo elementLooping3
enddo elementLooping3 !$OMP END PARALLEL DO
!$OMP END PARALLEL DO
else
print'(/,a,/)', ' << HOMOG >> Material Point terminally ill'
endif
end subroutine materialpoint_stressAndItsTangent2 end subroutine materialpoint_stressAndItsTangent2

View File

@ -163,8 +163,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
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]) if (.not. terminallyIll) &
cutBack = .false. ! reset cutBack status call materialpoint_stressAndItsTangent2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])
cutBack = .false.
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)