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

continue with cell based map (solver + homog)

This commit is contained in:
Sharan Roongta 2021-02-23 15:10:34 +01:00
parent d8112cc2e1
commit 9679c5403e
7 changed files with 73 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/grid/grid_damage_spectral.f90
View File

@ -160,7 +160,7 @@ function grid_damage_spectral_solution(timeinc) result(solution)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc !< increment in time for current solution timeinc !< increment in time for current solution
integer :: i, j, k, cell integer :: i, j, k, ce
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
PetscReal :: phi_min, phi_max, stagNorm, solnNorm PetscReal :: phi_min, phi_max, stagNorm, solnNorm
@ -193,10 +193,10 @@ function grid_damage_spectral_solution(timeinc) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updating damage state ! updating damage state
cell = 0 ce = 0
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)
cell = cell + 1 ce = ce + 1
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce)
enddo; enddo; enddo enddo; enddo; enddo
call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
@ -216,7 +216,7 @@ end function grid_damage_spectral_solution
subroutine grid_damage_spectral_forward(cutBack) subroutine grid_damage_spectral_forward(cutBack)
logical, intent(in) :: cutBack logical, intent(in) :: cutBack
integer :: i, j, k, cell integer :: i, j, k, ce
DM :: dm_local DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -226,14 +226,14 @@ subroutine grid_damage_spectral_forward(cutBack)
phi_stagInc = phi_lastInc phi_stagInc = phi_lastInc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reverting damage field state ! reverting damage field state
cell = 0 ce = 0
call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
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)
cell = cell + 1 ce = ce + 1
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce)
enddo; enddo; enddo enddo; enddo; enddo
else else
phi_lastInc = phi_current phi_lastInc = phi_current
@ -258,7 +258,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
f_scal f_scal
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: i, j, k, cell integer :: i, j, k, ce
real(pReal) :: phiDot, dPhiDot_dPhi, mobility real(pReal) :: phiDot, dPhiDot_dPhi, mobility
phi_current = x_scal phi_current = x_scal
@ -269,20 +269,20 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
call utilities_FFTscalarForward call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of damage field call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward call utilities_FFTvectorBackward
cell = 0 ce = 0
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)
cell = cell + 1 ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(1,cell) - K_ref, & vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(ce) - K_ref, &
vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward call utilities_FFTscalarBackward
cell = 0 ce = 0
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)
cell = cell + 1 ce = ce + 1
call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k), 1, cell) call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k),ce)
mobility = damage_nonlocal_getMobility(1,cell) mobility = damage_nonlocal_getMobility(ce)
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) & scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) &
+ mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
+ mu_ref*phi_current(i,j,k) + mu_ref*phi_current(i,j,k)
@ -310,15 +310,15 @@ end subroutine formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine updateReference subroutine updateReference
integer :: i,j,k,cell,ierr integer :: i,j,k,ce,ierr
cell = 0 ce = 0
K_ref = 0.0_pReal K_ref = 0.0_pReal
mu_ref = 0.0_pReal mu_ref = 0.0_pReal
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)
cell = cell + 1 ce = ce + 1
K_ref = K_ref + damage_nonlocal_getDiffusion(1,cell) K_ref = K_ref + damage_nonlocal_getDiffusion(ce)
mu_ref = mu_ref + damage_nonlocal_getMobility(1,cell) mu_ref = mu_ref + damage_nonlocal_getMobility(ce)
enddo; enddo; enddo enddo; enddo; enddo
K_ref = K_ref*wgt K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)

4
src/grid/grid_thermal_spectral.f90
View File

@ -268,7 +268,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = 0 ce = 0
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)
ce = ce + 1 ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(1,ce) - K_ref, & vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(ce) - K_ref, &
vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward call utilities_FFTvectorForward
@ -310,7 +310,7 @@ subroutine updateReference
mu_ref = 0.0_pReal mu_ref = 0.0_pReal
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)
ce = ce + 1 ce = ce + 1
K_ref = K_ref + thermal_conduction_getConductivity(1,ce) K_ref = K_ref + thermal_conduction_getConductivity(ce)
mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce) mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce)
enddo; enddo; enddo enddo; enddo; enddo
K_ref = K_ref*wgt K_ref = K_ref*wgt

47
src/homogenization.f90
View File

@ -133,11 +133,9 @@ module homogenization
end function mechanical_updateState end function mechanical_updateState
module function thermal_conduction_getConductivity(ip,el) result(K) module function thermal_conduction_getConductivity(ce) result(K)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: K real(pReal), dimension(3,3) :: K
end function thermal_conduction_getConductivity end function thermal_conduction_getConductivity
@ -169,23 +167,19 @@ module homogenization
module subroutine thermal_conduction_getSource(Tdot, ip, el) module subroutine thermal_conduction_getSource(Tdot, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, &
el !< element number el
real(pReal), intent(out) :: Tdot real(pReal), intent(out) :: Tdot
end subroutine thermal_conduction_getSource end subroutine thermal_conduction_getSource
module function damage_nonlocal_getMobility(ip,el) result(M) module function damage_nonlocal_getMobility(ce) result(M)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal) :: M real(pReal) :: M
end function damage_nonlocal_getMobility end function damage_nonlocal_getMobility
module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
real(pReal) :: & real(pReal) :: &
@ -193,11 +187,9 @@ module homogenization
end subroutine damage_nonlocal_getSourceAndItsTangent end subroutine damage_nonlocal_getSourceAndItsTangent
module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
@ -521,28 +513,25 @@ end subroutine damage_nonlocal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized non local damage diffusion tensor in reference configuration !> @brief returns homogenized non local damage diffusion tensor in reference configuration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function damage_nonlocal_getDiffusion(ip,el) function damage_nonlocal_getDiffusion(ce)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
damage_nonlocal_getDiffusion damage_nonlocal_getDiffusion
integer :: & integer :: &
homog, & ho, &
grain, & grain
ce
homog = material_homogenizationAt(el) ho = material_homogenizationAt2(ce)
damage_nonlocal_getDiffusion = 0.0_pReal damage_nonlocal_getDiffusion = 0.0_pReal
ce = (el-1)*discretization_nIPs + ip
do grain = 1, homogenization_Nconstituents(homog) do grain = 1, homogenization_Nconstituents(ho)
damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(grain,ce,lattice_D(1:3,1:3,material_phaseAt2(grain,ce))) crystallite_push33ToRef(grain,ce,lattice_D(1:3,1:3,material_phaseAt2(grain,ce)))
enddo enddo
damage_nonlocal_getDiffusion = & damage_nonlocal_getDiffusion = &
num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal) num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(ho),pReal)
end function damage_nonlocal_getDiffusion end function damage_nonlocal_getDiffusion

40
src/homogenization_damage.f90
View File

@ -85,22 +85,20 @@ end subroutine damage_partition
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Returns homogenized nonlocal damage mobility !> @brief Returns homogenized nonlocal damage mobility
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function damage_nonlocal_getMobility(ip,el) result(M) module function damage_nonlocal_getMobility(ce) result(M)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
integer :: & integer :: &
co co
real(pReal) :: M real(pReal) :: M
M = 0.0_pReal M = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
M = M + lattice_M(material_phaseAt(co,el)) M = M + lattice_M(material_phaseAt2(co,ce))
enddo enddo
M = M/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) M = M/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
end function damage_nonlocal_getMobility end function damage_nonlocal_getMobility
@ -108,11 +106,9 @@ end function damage_nonlocal_getMobility
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates homogenized damage driving forces !> @brief calculates homogenized damage driving forces
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
real(pReal) :: & real(pReal) :: &
@ -121,9 +117,9 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent end subroutine damage_nonlocal_getSourceAndItsTangent
@ -131,20 +127,18 @@ end subroutine damage_nonlocal_getSourceAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief updated nonlocal damage field with solution from damage phase field PDE !> @brief updated nonlocal damage field with solution from damage phase field PDE
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer :: & integer :: &
homog, & ho, &
offset me
homog = material_homogenizationAt(el) ho = material_homogenizationAt2(ce)
offset = material_homogenizationMemberAt(ip,el) me = material_homogenizationMemberAt2(ce)
damagestate_h(homog)%state(1,offset) = phi damagestate_h(ho)%state(1,me) = phi
end subroutine damage_nonlocal_putNonLocalDamage end subroutine damage_nonlocal_putNonLocalDamage

14
src/homogenization_thermal.f90
View File

@ -99,20 +99,16 @@ end subroutine thermal_homogenize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return homogenized thermal conductivity in reference configuration !> @brief return homogenized thermal conductivity in reference configuration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function thermal_conduction_getConductivity(ip,el) result(K) module function thermal_conduction_getConductivity(ce) result(K)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: K real(pReal), dimension(3,3) :: K
integer :: & integer :: &
co, & co
ce
K = 0.0_pReal K = 0.0_pReal
ce = (el-1)*discretization_nIPs + ip
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce))) K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce)))
enddo enddo
@ -223,8 +219,8 @@ end function homogenization_thermal_T
module subroutine thermal_conduction_getSource(Tdot, ip, el) module subroutine thermal_conduction_getSource(Tdot, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, &
el !< element number el
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
Tdot Tdot

6
src/phase.f90
View File

@ -231,10 +231,8 @@ module phase
end function phase_homogenizedC end function phase_homogenizedC
module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi !< damage parameter phi !< damage parameter
real(pReal), intent(inout) :: & real(pReal), intent(inout) :: &

12
src/phase_damage.f90
View File

@ -179,11 +179,9 @@ end subroutine damage_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief returns local part of nonlocal damage driving force !< @brief returns local part of nonlocal damage driving force
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi !< damage parameter phi !< damage parameter
real(pReal), intent(inout) :: & real(pReal), intent(inout) :: &
@ -201,9 +199,9 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi,
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
ph = material_phaseAt(co,el) ph = material_phaseAt2(co,ce)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt2(co,ce)
select case(phase_source(ph)) select case(phase_source(ph))
case (DAMAGE_ISOBRITTLE_ID) case (DAMAGE_ISOBRITTLE_ID)