(incomplete) split: spectral/FFT vs mech functionality

This commit is contained in:
Martin Diehl 2023-12-24 07:31:02 +01:00
parent 395b6b5438
commit 0d747ae5e4
No known key found for this signature in database
GPG Key ID: 1FD50837275A0A9B
8 changed files with 269 additions and 233 deletions

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -75,19 +75,6 @@ module spectral_utilities
termIll = .false. termIll = .false.
end type tSolutionState 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 type :: tNumerics
integer :: & integer :: &
divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints] 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_curlRMS, &
utilities_scalarGradient, & utilities_scalarGradient, &
utilities_vectorDivergence, & utilities_vectorDivergence, &
utilities_maskedCompliance, &
utilities_constitutiveResponse, &
utilities_calculateRate, &
utilities_forwardTensorField, &
utilities_updateCoords utilities_updateCoords
contains contains
@ -653,65 +636,6 @@ real(pREAL) function utilities_curlRMS(tensorField)
end function utilities_curlRMS 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. !> @brief Calculate gradient of scalar field.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -755,147 +679,6 @@ function utilities_vectorDivergence(field) result(div)
end function utilities_vectorDivergence 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. !> @brief Calculate Filter for Fourier convolution.
!> @details this is the full operator to calculate derivatives, i.e. 2 \pi i k for the !> @details this is the full operator to calculate derivatives, i.e. 2 \pi i k for the