diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index bce0e6db6..c8db7d8df 100644 --- a/src/Marc/materialpoint_Marc.f90 +++ b/src/Marc/materialpoint_Marc.f90 @@ -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) & - * homogenization_dPdF(i,m,k,n,ce) & - - math_delta(j,l) * homogenization_F(i,m,ce) * homogenization_P(k,m,ce) & + + materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) & + * homogenization_dPdF(i,m,k,n,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 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index cd6a60abe..3fd220ce5 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -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 diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 648f9ebbb..cce869653 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -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 diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index f7346055f..59835f250 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -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 diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 3fad9d59c..0e8ba4841 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -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 diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index 3b7cbda14..b5cc0b967 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -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 diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 new file mode 100644 index 000000000..0ea269319 --- /dev/null +++ b/src/grid/grid_mech_utilities.f90 @@ -0,0 +1,253 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief Utilities used by the mech grid solver variants +!-------------------------------------------------------------------------------------------------- +module grid_mech_utilities +#include + 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 !< - 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 diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index c178c5810..c8e638207 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -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 diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 19f497d57..4ea53d038 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -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 !< - 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