Merge branch 'solver-cleanup' into 'development'

Solver cleanup

See merge request damask/DAMASK!886
This commit is contained in:
Martin Diehl 2023-12-28 18:34:15 +00:00
commit f89e33a726
9 changed files with 279 additions and 240 deletions

View File

@ -30,7 +30,8 @@ module materialpoint_Marc
real(pREAL), dimension (:,:,:), allocatable, private :: &
materialpoint_cs !< Cauchy stress
real(pREAL), dimension (:,:,:,:), allocatable, private :: &
materialpoint_dcsdE !< Cauchy stress tangent
materialpoint_dcsdE, & !< Cauchy stress tangent
materialpoint_F !< deformation gradient
real(pREAL), dimension (:,:,:,:), allocatable, private :: &
materialpoint_dcsdE_knownGood !< known good tangent
@ -95,6 +96,7 @@ subroutine materialpoint_init()
print'(/,1x,a)', '<<<+- materialpoint init -+>>>'; flush(IO_STDOUT)
allocate(materialpoint_F( 3,3,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL)
allocate(materialpoint_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL)
allocate(materialpoint_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL)
allocate(materialpoint_dcsdE_knownGood(6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL)
@ -140,9 +142,10 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
if (iand(mode, materialpoint_RESTOREJACOBIAN) /= 0) &
materialpoint_dcsde = materialpoint_dcsde_knownGood
if (iand(mode, materialpoint_AGERESULTS) /= 0) call materialpoint_forward
if (iand(mode, materialpoint_AGERESULTS) /= 0) call materialpoint_forward()
homogenization_F(1:3,1:3,ce) = ffn1
materialpoint_F(1:3,1:3,ip,elCP) = ffn1
if (iand(mode, materialpoint_CALCRESULTS) /= 0) then
@ -167,17 +170,17 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
else terminalIllness
! translate from P to sigma
Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(homogenization_F(1:3,1:3,ce)))
J_inverse = 1.0_pREAL / math_det33(homogenization_F(1:3,1:3,ce))
Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pREAL / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
materialpoint_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
! translate from dP/dF to dCS/dE
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
H(i,j,k,l) = H(i,j,k,l) &
+ homogenization_F(j,m,ce) * homogenization_F(l,n,ce) &
+ materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) &
* homogenization_dPdF(i,m,k,n,ce) &
- math_delta(j,l) * homogenization_F(i,m,ce) * homogenization_P(k,m,ce) &
- math_delta(j,l) * materialpoint_F(i,m,ip,elCP) * homogenization_P(k,m,ce) &
+ 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))
end do; end do; end do; end do; end do; end do

View File

@ -23,6 +23,7 @@ program DAMASK_grid
use materialpoint
use material
use spectral_utilities
use grid_mech_utilities
use grid_mechanical_spectral_basic
use grid_mechanical_spectral_polarization
use grid_mechanical_FEM

View File

@ -44,7 +44,6 @@ module grid_damage_spectral
type(tNumerics) :: num
type(tSolutionParams) :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
SNES :: SNES_damage
@ -57,7 +56,7 @@ module grid_damage_spectral
! reference diffusion tensor, mobility etc.
integer :: totalIter = 0 !< total iteration in current increment
real(pREAL), dimension(3,3) :: K_ref
real(pREAL) :: mu_ref
real(pREAL) :: mu_ref, Delta_t_
public :: &
grid_damage_spectral_init, &
@ -207,8 +206,7 @@ end subroutine grid_damage_spectral_init
!--------------------------------------------------------------------------------------------------
function grid_damage_spectral_solution(Delta_t) result(solution)
real(pREAL), intent(in) :: &
Delta_t !< increment in time for current solution
real(pREAL), intent(in) :: Delta_t !< increment in time for current solution
type(tSolutionState) :: solution
PetscInt :: devNull
@ -222,7 +220,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%Delta_t = Delta_t
Delta_t_ = Delta_t
call SNESSolve(SNES_damage,PETSC_NULL_VEC,phi_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
@ -350,12 +348,12 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc)
ce = 0
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_phi(phi(i,j,k),ce)) &
r(i,j,k) = Delta_t_*(r(i,j,k) + homogenization_f_phi(phi(i,j,k),ce)) &
+ homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi(i,j,k)) &
+ mu_ref*phi(i,j,k)
end do; end do; end do
r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t),phi_lastInc),num%phi_min) &
r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, Delta_t_),phi_lastInc),num%phi_min) &
- phi
end associate
err_PETSc = 0

View File

@ -23,6 +23,7 @@ module grid_mechanical_FEM
use math
use rotations
use spectral_utilities
use grid_mech_utilities
use config
use homogenization
use discretization

View File

@ -23,6 +23,7 @@ module grid_mechanical_spectral_basic
use math
use rotations
use spectral_utilities
use grid_mech_utilities
use homogenization
use discretization_grid

View File

@ -23,6 +23,7 @@ module grid_mechanical_spectral_polarization
use math
use rotations
use spectral_utilities
use grid_mech_utilities
use config
use homogenization
use discretization_grid

View File

@ -0,0 +1,253 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Utilities used by the mech grid solver variants
!--------------------------------------------------------------------------------------------------
module grid_mech_utilities
#include <petsc/finclude/petscsys.h>
use PETScSys
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
use MPI_f08
#endif
use prec
use parallelization
use math
use rotations
use IO
use discretization_grid
use discretization
use spectral_utilities
use homogenization
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
implicit none(type,external)
#else
implicit none
#endif
private
!--------------------------------------------------------------------------------------------------
! derived types
type, public :: tBoundaryCondition !< set of parameters defining a boundary condition
real(pREAL), dimension(3,3) :: values = 0.0_pREAL
logical, dimension(3,3) :: mask = .true.
character(len=:), allocatable :: myType
end type tBoundaryCondition
type, public :: tSolutionParams
real(pREAL), dimension(3,3) :: stress_BC
logical, dimension(3,3) :: stress_mask
type(tRotation) :: rotation_BC
real(pREAL) :: Delta_t
end type tSolutionParams
public :: &
utilities_maskedCompliance, &
utilities_constitutiveResponse, &
utilities_calculateRate, &
utilities_forwardTensorField
contains
!--------------------------------------------------------------------------------------------------
!> @brief Calculate masked compliance tensor used to adjust F to fullfill stress BC.
!--------------------------------------------------------------------------------------------------
function utilities_maskedCompliance(rot_BC,mask_stress,C)
real(pREAL), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance
real(pREAL), intent(in), dimension(3,3,3,3) :: C !< current average stiffness
type(tRotation), intent(in) :: rot_BC !< rotation of load frame
logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC
integer :: i, j
logical, dimension(9) :: mask_stressVector
logical, dimension(9,9) :: mask
real(pREAL), dimension(9,9) :: temp99_real
integer :: size_reduced = 0
real(pREAL), dimension(:,:), allocatable :: &
s_reduced, & !< reduced compliance matrix (depending on number of stress BC)
c_reduced, & !< reduced stiffness (depending on number of stress BC)
sTimesC !< temp variable to check inversion
logical :: errmatinv
character(len=pSTRLEN):: formatString
mask_stressVector = .not. reshape(transpose(mask_stress), [9])
size_reduced = count(mask_stressVector)
if (size_reduced > 0) then
temp99_real = math_3333to99(rot_BC%rotate(C))
do i = 1,9; do j = 1,9
mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j)
end do; end do
c_reduced = reshape(pack(temp99_Real,mask),[size_reduced,size_reduced])
allocate(s_reduced,mold = c_reduced)
call math_invert(s_reduced, errmatinv, c_reduced) ! invert reduced stiffness
if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
!--------------------------------------------------------------------------------------------------
! check if inversion was successful
sTimesC = matmul(c_reduced,s_reduced)
errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pREAL))
if (errmatinv) then
write(formatString, '(i2)') size_reduced
formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced))
print trim(formatString), 'C (load) ', transpose(c_reduced)
print trim(formatString), 'S (load) ', transpose(s_reduced)
if (errmatinv) error stop 'matrix inversion error'
end if
temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pREAL),[9,9])
else
temp99_real = 0.0_pREAL
end if
utilities_maskedCompliance = math_99to3333(temp99_Real)
end function utilities_maskedCompliance
!--------------------------------------------------------------------------------------------------
!> @brief Calculate constitutive response.
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
F,Delta_t,rotation_BC)
real(pREAL), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pREAL), intent(out), dimension(3,3) :: P_av !< average PK stress
real(pREAL), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress
real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: F !< deformation gradient target
real(pREAL), intent(in) :: Delta_t !< loading time
type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame
integer :: i
integer(MPI_INTEGER_KIND) :: err_MPI
real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min
real(pREAL) :: dPdF_norm_max, dPdF_norm_min
real(pREAL), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF
print'(/,1x,a)', '... evaluating constitutive response ......................................'
flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field
call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3) ! calculate P field
if (.not. terminallyIll) &
call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3)
if (.not. terminallyIll) &
call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (present(rotation_BC)) then
if (any(dNeq(rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pREAL)))) &
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
'Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pREAL
P_av = rotation_BC%rotate(P_av)
end if
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
'Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pREAL
flush(IO_STDOUT)
dPdF_max = 0.0_pREAL
dPdF_norm_max = 0.0_pREAL
dPdF_min = huge(1.0_pREAL)
dPdF_norm_min = huge(1.0_pREAL)
do i = 1, product(cells(1:2))*cells3
if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
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,i)**2)
end if
if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
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,i)**2)
end if
end do
valueAndRank = [dPdF_norm_max,real(worldrank,pREAL)]
call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MAXLOC,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Bcast(dPdF_max,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
valueAndRank = [dPdF_norm_min,real(worldrank,pREAL)]
call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MINLOC,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Bcast(dPdF_min,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
C_minmaxAvg = 0.5_pREAL*(dPdF_max + dPdF_min)
C_volAvg = sum(homogenization_dPdF,dim=5)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
C_volAvg = C_volAvg * wgt
end subroutine utilities_constitutiveResponse
!--------------------------------------------------------------------------------------------------
!> @brief Calculate forward rate, either as local guess or as homogeneous add on.
!--------------------------------------------------------------------------------------------------
pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
real(pREAL), intent(in), dimension(3,3) :: &
avRate !< homogeneous addon
real(pREAL), intent(in) :: &
dt !< Delta_t between field0 and field
logical, intent(in) :: &
heterogeneous !< calculate field of rates
real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
field0, & !< data of previous step
field !< data of current step
real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_calculateRate
utilities_calculateRate = merge((field-field0) / dt, &
spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3), &
heterogeneous)
end function utilities_calculateRate
!--------------------------------------------------------------------------------------------------
!> @brief forwards a field with a pointwise given rate, if aim is given,
!> ensures that the average matches the aim
!--------------------------------------------------------------------------------------------------
function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim)
real(pREAL), intent(in) :: &
Delta_t !< Delta_t of current step
real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
field_lastInc, & !< initial field
rate !< rate by which to forward
real(pREAL), intent(in), optional, dimension(3,3) :: &
aim !< average field value aim
real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_forwardTensorField
real(pREAL), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
integer(MPI_INTEGER_KIND) :: err_MPI
utilities_forwardTensorField = field_lastInc + rate*Delta_t
if (present(aim)) then !< correct to match average
fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt
call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
fieldDiff = fieldDiff - aim
utilities_forwardTensorField = utilities_forwardTensorField &
- spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
end if
end function utilities_forwardTensorField
end module grid_mech_utilities

View File

@ -43,7 +43,6 @@ module grid_thermal_spectral
type(tNumerics) :: num
type(tSolutionParams) :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
SNES :: SNES_thermal
@ -56,7 +55,7 @@ module grid_thermal_spectral
! reference diffusion tensor, mobility etc.
integer :: totalIter = 0 !< total iteration in current increment
real(pREAL), dimension(3,3) :: K_ref
real(pREAL) :: mu_ref
real(pREAL) :: mu_ref, Delta_t_
public :: &
grid_thermal_spectral_init, &
@ -186,8 +185,7 @@ end subroutine grid_thermal_spectral_init
!--------------------------------------------------------------------------------------------------
function grid_thermal_spectral_solution(Delta_t) result(solution)
real(pREAL), intent(in) :: &
Delta_t !< increment in time for current solution
real(pREAL), intent(in) :: Delta_t !< increment in time for current solution
type(tSolutionState) :: solution
PetscInt :: devNull
@ -201,7 +199,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%Delta_t = Delta_t
Delta_t_ = Delta_t
call SNESSolve(SNES_thermal,PETSC_NULL_VEC,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
@ -227,7 +225,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
T_stagInc = T
call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), &
reshape(T-T_lastInc,[product(cells(1:2))*cells3])/params%Delta_t)
reshape(T-T_lastInc,[product(cells(1:2))*cells3])/Delta_t_)
call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc)
CHKERRQ(err_PETSc)
@ -264,7 +262,7 @@ subroutine grid_thermal_spectral_forward(cutBack)
T = T_lastInc
T_stagInc = T_lastInc
else
dotT_lastInc = (T - T_lastInc)/params%Delta_t
dotT_lastInc = (T - T_lastInc)/Delta_t_
T_lastInc = T
call updateReference()
end if
@ -336,13 +334,13 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc)
ce = 0
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_T(ce)) &
r(i,j,k) = Delta_t_*(r(i,j,k) + homogenization_f_T(ce)) &
+ homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T(i,j,k)) &
+ mu_ref*T(i,j,k)
end do; end do; end do
r = T &
- utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t)
- utilities_GreenConvolution(r, K_ref, mu_ref, Delta_t_)
end associate
err_PETSc = 0

View File

@ -75,19 +75,6 @@ module spectral_utilities
termIll = .false.
end type tSolutionState
type, public :: tBoundaryCondition !< set of parameters defining a boundary condition
real(pREAL), dimension(3,3) :: values = 0.0_pREAL
logical, dimension(3,3) :: mask = .true.
character(len=:), allocatable :: myType
end type tBoundaryCondition
type, public :: tSolutionParams
real(pREAL), dimension(3,3) :: stress_BC
logical, dimension(3,3) :: stress_mask
type(tRotation) :: rotation_BC
real(pREAL) :: Delta_t
end type tSolutionParams
type :: tNumerics
integer :: &
divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints]
@ -121,10 +108,6 @@ module spectral_utilities
utilities_curlRMS, &
utilities_scalarGradient, &
utilities_vectorDivergence, &
utilities_maskedCompliance, &
utilities_constitutiveResponse, &
utilities_calculateRate, &
utilities_forwardTensorField, &
utilities_updateCoords
contains
@ -653,65 +636,6 @@ real(pREAL) function utilities_curlRMS(tensorField)
end function utilities_curlRMS
!--------------------------------------------------------------------------------------------------
!> @brief Calculate masked compliance tensor used to adjust F to fullfill stress BC.
!--------------------------------------------------------------------------------------------------
function utilities_maskedCompliance(rot_BC,mask_stress,C)
real(pREAL), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance
real(pREAL), intent(in), dimension(3,3,3,3) :: C !< current average stiffness
type(tRotation), intent(in) :: rot_BC !< rotation of load frame
logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC
integer :: i, j
logical, dimension(9) :: mask_stressVector
logical, dimension(9,9) :: mask
real(pREAL), dimension(9,9) :: temp99_real
integer :: size_reduced = 0
real(pREAL), dimension(:,:), allocatable :: &
s_reduced, & !< reduced compliance matrix (depending on number of stress BC)
c_reduced, & !< reduced stiffness (depending on number of stress BC)
sTimesC !< temp variable to check inversion
logical :: errmatinv
character(len=pSTRLEN):: formatString
mask_stressVector = .not. reshape(transpose(mask_stress), [9])
size_reduced = count(mask_stressVector)
if (size_reduced > 0) then
temp99_real = math_3333to99(rot_BC%rotate(C))
do i = 1,9; do j = 1,9
mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j)
end do; end do
c_reduced = reshape(pack(temp99_Real,mask),[size_reduced,size_reduced])
allocate(s_reduced,mold = c_reduced)
call math_invert(s_reduced, errmatinv, c_reduced) ! invert reduced stiffness
if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
!--------------------------------------------------------------------------------------------------
! check if inversion was successful
sTimesC = matmul(c_reduced,s_reduced)
errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pREAL))
if (errmatinv) then
write(formatString, '(i2)') size_reduced
formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced))
print trim(formatString), 'C (load) ', transpose(c_reduced)
print trim(formatString), 'S (load) ', transpose(s_reduced)
if (errmatinv) error stop 'matrix inversion error'
end if
temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pREAL),[9,9])
else
temp99_real = 0.0_pREAL
end if
utilities_maskedCompliance = math_99to3333(temp99_Real)
end function utilities_maskedCompliance
!--------------------------------------------------------------------------------------------------
!> @brief Calculate gradient of scalar field.
!--------------------------------------------------------------------------------------------------
@ -755,147 +679,6 @@ function utilities_vectorDivergence(field) result(div)
end function utilities_vectorDivergence
!--------------------------------------------------------------------------------------------------
!> @brief Calculate constitutive response.
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
F,Delta_t,rotation_BC)
real(pREAL), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pREAL), intent(out), dimension(3,3) :: P_av !< average PK stress
real(pREAL), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress
real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: F !< deformation gradient target
real(pREAL), intent(in) :: Delta_t !< loading time
type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame
integer :: i
integer(MPI_INTEGER_KIND) :: err_MPI
real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min
real(pREAL) :: dPdF_norm_max, dPdF_norm_min
real(pREAL), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF
print'(/,1x,a)', '... evaluating constitutive response ......................................'
flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field
call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3) ! calculate P field
if (.not. terminallyIll) &
call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3)
if (.not. terminallyIll) &
call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (present(rotation_BC)) then
if (any(dNeq(rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pREAL)))) &
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
'Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pREAL
P_av = rotation_BC%rotate(P_av)
end if
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
'Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pREAL
flush(IO_STDOUT)
dPdF_max = 0.0_pREAL
dPdF_norm_max = 0.0_pREAL
dPdF_min = huge(1.0_pREAL)
dPdF_norm_min = huge(1.0_pREAL)
do i = 1, product(cells(1:2))*cells3
if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
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,i)**2)
end if
if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
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,i)**2)
end if
end do
valueAndRank = [dPdF_norm_max,real(worldrank,pREAL)]
call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MAXLOC,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Bcast(dPdF_max,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
valueAndRank = [dPdF_norm_min,real(worldrank,pREAL)]
call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MINLOC,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call MPI_Bcast(dPdF_min,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
C_minmaxAvg = 0.5_pREAL*(dPdF_max + dPdF_min)
C_volAvg = sum(homogenization_dPdF,dim=5)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
C_volAvg = C_volAvg * wgt
end subroutine utilities_constitutiveResponse
!--------------------------------------------------------------------------------------------------
!> @brief Calculate forward rate, either as local guess or as homogeneous add on.
!--------------------------------------------------------------------------------------------------
pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
real(pREAL), intent(in), dimension(3,3) :: &
avRate !< homogeneous addon
real(pREAL), intent(in) :: &
dt !< Delta_t between field0 and field
logical, intent(in) :: &
heterogeneous !< calculate field of rates
real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
field0, & !< data of previous step
field !< data of current step
real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_calculateRate
utilities_calculateRate = merge((field-field0) / dt, &
spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3), &
heterogeneous)
end function utilities_calculateRate
!--------------------------------------------------------------------------------------------------
!> @brief forwards a field with a pointwise given rate, if aim is given,
!> ensures that the average matches the aim
!--------------------------------------------------------------------------------------------------
function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim)
real(pREAL), intent(in) :: &
Delta_t !< Delta_t of current step
real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
field_lastInc, & !< initial field
rate !< rate by which to forward
real(pREAL), intent(in), optional, dimension(3,3) :: &
aim !< average field value aim
real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_forwardTensorField
real(pREAL), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
integer(MPI_INTEGER_KIND) :: err_MPI
utilities_forwardTensorField = field_lastInc + rate*Delta_t
if (present(aim)) then !< correct to match average
fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt
call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
fieldDiff = fieldDiff - aim
utilities_forwardTensorField = utilities_forwardTensorField &
- spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
end if
end function utilities_forwardTensorField
!--------------------------------------------------------------------------------------------------
!> @brief Calculate Filter for Fourier convolution.
!> @details this is the full operator to calculate derivatives, i.e. 2 \pi i k for the