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

No FEsolving

This commit is contained in:
Martin Diehl 2021-01-06 14:24:02 +01:00 committed by Sharan Roongta
parent f8dd5df0cc
commit 437d91495b
17 changed files with 72 additions and 130 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
src/CPFEM.f90
View File

@ -5,7 +5,6 @@
!--------------------------------------------------------------------------------------------------
module CPFEM
use prec
use FEsolving
use math
use rotations
use YAML_types
@ -197,11 +196,9 @@ 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
FEsolving_execElem = elCP
FEsolving_execIP = ip
if (debugCPFEM%extensive) &
print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
call materialpoint_stressAndItsTangent(dt)
call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP])
terminalIllness: if (terminallyIll) then

1
src/CPFEM2.f90
View File

@ -6,7 +6,6 @@
module CPFEM2
use prec
use config
use FEsolving
use math
use rotations
use YAML_types

1
src/DAMASK_marc.f90
View File

@ -176,7 +176,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
use DAMASK_interface
use config
use YAML_types
use FEsolving
use discretization_marc
use homogenization
use CPFEM

15
src/FEsolving.f90
View File

@ -1,15 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief global variables for flow control
!--------------------------------------------------------------------------------------------------
module FEsolving
implicit none
public
integer, dimension(2) :: &
FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element
FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP
end module FEsolving

1
src/commercialFEM_fileList.f90
View File

@ -13,7 +13,6 @@
#include "math.f90"
#include "quaternions.f90"
#include "rotations.f90"
#include "FEsolving.f90"
#include "element.f90"
#include "HDF5_utilities.f90"
#include "results.f90"

55
src/constitutive.f90
View File

@ -16,7 +16,6 @@ module constitutive
use parallelization
use HDF5_utilities
use DAMASK_interface
use FEsolving
use results
implicit none
@ -940,8 +939,8 @@ subroutine crystallite_init
flush(IO_STDOUT)
!$OMP PARALLEL DO PRIVATE(ph,me)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1), FEsolving_execIP(2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
do el = 1, size(material_phaseMemberAt,3)
do ip = 1, size(material_phaseMemberAt,2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
@ -967,14 +966,14 @@ subroutine crystallite_init
crystallite_partitionedF0 = crystallite_F0
crystallite_partitionedF = crystallite_F0
call crystallite_orientations()
!$OMP PARALLEL DO PRIVATE(ph,me)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
do el = 1, size(material_phaseMemberAt,3)
do ip = 1, size(material_phaseMemberAt,2)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call crystallite_orientations(co,ip,el)
call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el) ! update dependent state variables to be consistent with basic states
enddo
enddo
@ -1089,7 +1088,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
el !< counter in element loop
integer :: &
o, &
p, pp, m
p, ph, me
real(pReal), dimension(3,3) :: devNull, &
invSubFp0,invSubFi0,invFp,invFi, &
@ -1109,19 +1108,19 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
real(pReal), dimension(9,9):: temp_99
logical :: error
pp = material_phaseAt(co,el)
m = material_phaseMemberAt(co,ip,el)
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
crystallite_Fe(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el)
call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, &
crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
co,ip,el)
invFp = math_inv33(constitutive_mech_Fp(pp)%data(1:3,1:3,m))
invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m))
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me))
invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el))
invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el))
@ -1150,7 +1149,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!--------------------------------------------------------------------------------------------------
@ -1210,34 +1209,20 @@ end function crystallite_stressTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates orientations
!--------------------------------------------------------------------------------------------------
subroutine crystallite_orientations
subroutine crystallite_orientations(co,ip,el)
integer &
integer, intent(in) :: &
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el !< counter in element loop
!$OMP PARALLEL DO
do el = FEsolving_execElem(1),FEsolving_execElem(2)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el))))
enddo; enddo; enddo
!$OMP END PARALLEL DO
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el))))
if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(1,el)),ip,el)
nonlocalPresent: if (any(plasticState%nonlocal)) then
!$OMP PARALLEL DO
do el = FEsolving_execElem(1),FEsolving_execElem(2)
if (plasticState(material_phaseAt(1,el))%nonlocal) then
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(1,el)),ip,el)
enddo
endif
enddo
!$OMP END PARALLEL DO
endif nonlocalPresent
end subroutine crystallite_orientations

4
src/grid/discretization_grid.f90
View File

@ -19,7 +19,6 @@ module discretization_grid
use results
use discretization
use geometry_plastic_nonlocal
use FEsolving
implicit none
private
@ -117,9 +116,6 @@ subroutine discretization_grid_init(restart)
(grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process
worldrank+1==worldsize))
FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements
FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP
!--------------------------------------------------------------------------------------------------
! store geometry information for post processing
if(.not. restart) then

1
src/grid/grid_mech_FEM.f90
View File

@ -18,7 +18,6 @@ module grid_mech_FEM
use math
use rotations
use spectral_utilities
use FEsolving
use config
use homogenization
use discretization

1
src/grid/grid_mech_spectral_basic.f90
View File

@ -18,7 +18,6 @@ module grid_mech_spectral_basic
use math
use rotations
use spectral_utilities
use FEsolving
use config
use homogenization
use discretization_grid

1
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -18,7 +18,6 @@ module grid_mech_spectral_polarisation
use math
use rotations
use spectral_utilities
use FEsolving
use config
use homogenization
use discretization_grid

4
src/grid/spectral_utilities.f90
View File

@ -810,9 +810,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
print'(/,a)', ' ... evaluating constitutive response ......................................'
flush(IO_STDOUT)
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) ! calculate P field
call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field
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

76
src/homogenization.f90
View File

@ -11,7 +11,6 @@ module homogenization
use math
use material
use constitutive
use FEsolving
use discretization
use thermal_isothermal
use thermal_conduction
@ -144,27 +143,29 @@ end subroutine homogenization_init
!--------------------------------------------------------------------------------------------------
!> @brief parallelized calculation of stress and corresponding tangent at material points
!--------------------------------------------------------------------------------------------------
subroutine materialpoint_stressAndItsTangent(dt)
subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem)
real(pReal), intent(in) :: dt !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer :: &
NiterationHomog, &
NiterationMPstate, &
ip, & !< integration point number
el, & !< element number
myNgrains, co, ce
myNgrains, co, ce, ho
real(pReal) :: &
subFrac, &
subStep
logical :: &
requested, &
converged
logical, dimension(2) :: &
doneAndHappy
!$OMP PARALLEL DO PRIVATE(ce,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,requested,doneAndHappy)
!$OMP PARALLEL DO PRIVATE(ce,ho,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
!--------------------------------------------------------------------------------------------------
@ -174,21 +175,19 @@ subroutine materialpoint_stressAndItsTangent(dt)
subFrac = 0.0_pReal
converged = .false. ! pretend failed step ...
subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
requested = .true. ! everybody requires calculation
if (homogState(material_homogenizationAt(el))%sizeState > 0) &
homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el))
if (homogState(ho)%sizeState > 0) &
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(ho)%State0( :,material_homogenizationMemberAt(ip,el))
if (damageState(ho)%sizeState > 0) &
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(ho)%State0( :,material_homogenizationMemberAt(ip,el))
if (damageState(material_homogenizationAt(el))%sizeState > 0) &
damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(material_homogenizationAt(el))%State0( :,material_homogenizationMemberAt(ip,el))
NiterationHomog = 0
cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
myNgrains = homogenization_Nconstituents(material_homogenizationAt(el))
if (converged) then
subFrac = subFrac + subStep
subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration
@ -198,22 +197,20 @@ subroutine materialpoint_stressAndItsTangent(dt)
! wind forward grain starting point
call constitutive_windForward(ip,el)
if(homogState(material_homogenizationAt(el))%sizeState > 0) &
homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el))
if(damageState(material_homogenizationAt(el))%sizeState > 0) &
damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(material_homogenizationAt(el))%State (:,material_homogenizationMemberAt(ip,el))
if(homogState(ho)%sizeState > 0) &
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(ho)%State (:,material_homogenizationMemberAt(ip,el))
if(damageState(ho)%sizeState > 0) &
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(ho)%State (:,material_homogenizationMemberAt(ip,el))
endif steppingNeeded
else
if ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense
if (.not. terminallyIll) then ! so first signals terminally ill...
if (.not. terminallyIll) & ! so first signals terminally ill...
print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
endif
terminallyIll = .true. ! ...and kills all others
else ! cutback makes sense
subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback
@ -221,23 +218,19 @@ subroutine materialpoint_stressAndItsTangent(dt)
call crystallite_restore(ip,el,subStep < 1.0_pReal)
call constitutive_restore(ip,el)
if(homogState(material_homogenizationAt(el))%sizeState > 0) &
homogState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = &
homogState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el))
if(damageState(material_homogenizationAt(el))%sizeState > 0) &
damageState(material_homogenizationAt(el))%State( :,material_homogenizationMemberAt(ip,el)) = &
damageState(material_homogenizationAt(el))%subState0(:,material_homogenizationMemberAt(ip,el))
if(homogState(ho)%sizeState > 0) &
homogState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = &
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el))
if(damageState(ho)%sizeState > 0) &
damageState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = &
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el))
endif
endif
if (subStep > num%subStepMinHomog) then
requested = .true.
doneAndHappy = [.false.,.true.]
endif
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
NiterationMPstate = 0
convergenceLooping: do while (.not. terminallyIll .and. requested &
convergenceLooping: do while (.not. terminallyIll &
.and. .not. doneAndHappy(1) &
.and. NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1
@ -245,7 +238,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
!--------------------------------------------------------------------------------------------------
! deformation partitioning
if(requested .and. .not. doneAndHappy(1)) then ! requested but not yet done
if (.not. doneAndHappy(1)) then
ce = (el-1)*discretization_nIPs + ip
call mech_partition(homogenization_F0(1:3,1:3,ce) &
+ (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))&
@ -255,10 +248,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
do co = 1, myNgrains
converged = converged .and. crystallite_stress(dt*subStep,co,ip,el)
enddo
endif
if (requested .and. .not. doneAndHappy(1)) then
if (.not. converged) then
doneAndHappy = [.true.,.false.]
else
@ -281,10 +271,14 @@ subroutine materialpoint_stressAndItsTangent(dt)
!$OMP END PARALLEL DO
if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(ho,myNgrains)
elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho)
IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
do co = 1, myNgrains
call crystallite_orientations(co,ip,el)
enddo
call mech_homogenize(ip,el)
enddo IpLooping3
enddo elementLooping3

24
src/homogenization_mech.f90
View File

@ -128,35 +128,35 @@ module subroutine mech_homogenize(ip,el)
integer, intent(in) :: &
ip, & !< integration point
el !< element number
integer :: c,m
integer :: co,ce
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
m = (el-1)* discretization_nIPs + ip
ce = (el-1)* discretization_nIPs + ip
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
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,m) = crystallite_stressTangent(1,ip,el)
homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
enddo
call mech_isostrain_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,m), &
homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
enddo
call mech_RGC_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,m), &
homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))

4
src/marc/discretization_marc.f90
View File

@ -12,7 +12,6 @@ module discretization_marc
use DAMASK_interface
use IO
use config
use FEsolving
use element
use discretization
use geometry_plastic_nonlocal
@ -89,9 +88,6 @@ subroutine discretization_marc_init
if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,ext_msg='element')
if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP')
FEsolving_execElem = [1,nElems]
FEsolving_execIP = [1,elem%nIPs]
allocate(cellNodeDefinition(elem%nNodes-1))
allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems))
call buildCells(connectivity_cell,cellNodeDefinition,&

1
src/mesh/DAMASK_mesh.f90
View File

@ -15,7 +15,6 @@ program DAMASK_mesh
use IO
use math
use CPFEM2
use FEsolving
use config
use discretization_mesh
use FEM_Utilities

2
src/mesh/FEM_utilities.f90
View File

@ -160,7 +160,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
print'(/,a)', ' ... evaluating constitutive response ......................................'
call materialpoint_stressAndItsTangent(timeinc) ! calculate P field
call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field
cutBack = .false. ! reset cutBack status

6
src/mesh/discretization_mesh.f90
View File

@ -18,7 +18,6 @@ module discretization_mesh
use config
use discretization
use results
use FEsolving
use FEM_quadrature
use YAML_types
use prec
@ -30,7 +29,7 @@ module discretization_mesh
mesh_Nboundaries, &
mesh_NcpElemsGlobal
integer :: &
integer, public, protected :: &
mesh_NcpElems !< total number of CP elements in mesh
!!!! BEGIN DEPRECATED !!!!!
@ -174,9 +173,6 @@ subroutine discretization_mesh_init(restart)
if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element')
if (debug_ip < 1 .or. debug_ip > mesh_maxNips) call IO_error(602,ext_msg='IP')
FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements
FEsolving_execIP = [1,mesh_maxNips]
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
call discretization_init(materialAt,&