avoid code duplication
This commit is contained in:
Martin Diehl
parent
83453d10ef
commit
8a9d3f8d6d
|
|
@ -465,7 +465,7 @@ program DAMASK_spectral
|
|||
cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
|
||||
deformation_BC = loadCases(currentLoadCase)%deformation, &
|
||||
stress_BC = loadCases(currentLoadCase)%stress, &
|
||||
rotation_BC = loadCases(currentLoadCase)%rot%asMatrix())
|
||||
rotation_BC = loadCases(currentLoadCase)%rot)
|
||||
|
||||
case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
|
||||
case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack)
|
||||
|
|
@ -484,7 +484,7 @@ program DAMASK_spectral
|
|||
solres(field) = mech_solution (&
|
||||
incInfo,timeinc,timeIncOld, &
|
||||
stress_BC = loadCases(currentLoadCase)%stress, &
|
||||
rotation_BC = loadCases(currentLoadCase)%rot%asMatrix())
|
||||
rotation_BC = loadCases(currentLoadCase)%rot)
|
||||
|
||||
case(FIELD_THERMAL_ID)
|
||||
solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld)
|
||||
|
|
|
|||
|
|
@ -242,7 +242,8 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
|
|||
timeinc_old !< time increment of last successful increment
|
||||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC
|
||||
real(pReal), dimension(3,3), intent(in) :: rotation_BC
|
||||
type(rotation), intent(in) :: &
|
||||
rotation_BC
|
||||
type(tSolutionState) :: &
|
||||
solution
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -254,7 +255,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! update stiffness (and gamma operator)
|
||||
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
|
||||
S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! set module wide available data
|
||||
params%stress_mask = stress_BC%maskFloat
|
||||
|
|
@ -297,7 +298,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
|
|||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC, &
|
||||
deformation_BC
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
type(rotation), intent(in) :: &
|
||||
rotation_BC
|
||||
PetscErrorCode :: ierr
|
||||
PetscScalar, pointer, dimension(:,:,:,:) :: &
|
||||
|
|
@ -482,7 +483,7 @@ subroutine formResidual(da_local,x_local, &
|
|||
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax
|
||||
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
||||
' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
|
||||
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
||||
' deformation gradient aim =', transpose(F_aim)
|
||||
flush(6)
|
||||
|
|
@ -498,7 +499,7 @@ subroutine formResidual(da_local,x_local, &
|
|||
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
|
||||
enddo; enddo; enddo
|
||||
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
|
||||
F(1:3,1:3,ii,jj,kk) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem))
|
||||
F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotTensor2(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
|
||||
enddo; enddo; enddo
|
||||
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
|
||||
|
||||
|
|
@ -506,7 +507,7 @@ subroutine formResidual(da_local,x_local, &
|
|||
! evaluate constitutive response
|
||||
call Utilities_constitutiveResponse(P_current,&
|
||||
P_av,C_volAvg,devNull, &
|
||||
F,params%timeinc,params%rotation_BC)
|
||||
F,params%timeinc,params%rotation_BC%asMatrix())
|
||||
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -212,7 +212,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
|
|||
timeinc_old !< time increment of last successful increment
|
||||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC
|
||||
real(pReal), dimension(3,3), intent(in) :: rotation_BC
|
||||
type(rotation), intent(in) :: &
|
||||
rotation_BC
|
||||
type(tSolutionState) :: &
|
||||
solution
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -224,7 +225,7 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! update stiffness (and gamma operator)
|
||||
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
|
||||
S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg)
|
||||
if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -269,7 +270,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
|
|||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC, &
|
||||
deformation_BC
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
type(rotation), intent(in) :: &
|
||||
rotation_BC
|
||||
PetscErrorCode :: ierr
|
||||
PetscScalar, dimension(:,:,:,:), pointer :: F
|
||||
|
|
@ -299,9 +300,9 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
|
|||
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
|
||||
endif
|
||||
|
||||
Fdot = utilities_calculateRate(guess, &
|
||||
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
|
||||
math_rotate_backward33(F_aimDot,rotation_BC))
|
||||
Fdot = utilities_calculateRate(guess, &
|
||||
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
|
||||
rotation_BC%rotTensor2(F_aimDot,active=.true.))
|
||||
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
|
||||
|
||||
materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
|
||||
|
|
@ -311,7 +312,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
|
|||
! update average and local deformation gradients
|
||||
F_aim = F_aim_lastInc + F_aimDot * timeinc
|
||||
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
|
||||
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
|
||||
rotation_BC%rotTensor2(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
|
||||
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
|
||||
|
||||
end subroutine grid_mech_spectral_basic_forward
|
||||
|
|
@ -446,7 +447,7 @@ subroutine formResidual(in, F, &
|
|||
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
|
||||
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
||||
' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
|
||||
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
||||
' deformation gradient aim =', transpose(F_aim)
|
||||
flush(6)
|
||||
|
|
@ -456,7 +457,7 @@ subroutine formResidual(in, F, &
|
|||
! evaluate constitutive response
|
||||
call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory)
|
||||
P_av,C_volAvg,C_minMaxAvg, &
|
||||
F,params%timeinc,params%rotation_BC)
|
||||
F,params%timeinc,params%rotation_BC%asMatrix())
|
||||
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -471,7 +472,7 @@ subroutine formResidual(in, F, &
|
|||
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform
|
||||
call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
|
||||
err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
|
||||
call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg
|
||||
call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier
|
||||
call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -14,6 +14,7 @@ module grid_mech_spectral_polarisation
|
|||
use DAMASK_interface
|
||||
use HDF5_utilities
|
||||
use math
|
||||
use rotations
|
||||
use spectral_utilities
|
||||
use IO
|
||||
use FEsolving
|
||||
|
|
@ -222,12 +223,13 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
|
|||
! input data for solution
|
||||
character(len=*), intent(in) :: &
|
||||
incInfoIn
|
||||
real(pReal), intent(in) :: &
|
||||
real(pReal), intent(in) :: &
|
||||
timeinc, & !< time increment of current solution
|
||||
timeinc_old !< time increment of last successful increment
|
||||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC
|
||||
real(pReal), dimension(3,3), intent(in) :: rotation_BC
|
||||
type(rotation), intent(in) :: &
|
||||
rotation_BC
|
||||
type(tSolutionState) :: &
|
||||
solution
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -239,7 +241,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! update stiffness (and gamma operator)
|
||||
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
|
||||
S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg)
|
||||
if (num%update_gamma) then
|
||||
call utilities_updateGamma(C_minMaxAvg)
|
||||
C_scale = C_minMaxAvg
|
||||
|
|
@ -288,7 +290,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
|
|||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC, &
|
||||
deformation_BC
|
||||
real(pReal), dimension(3,3), intent(in) ::&
|
||||
type(rotation), intent(in) :: &
|
||||
rotation_BC
|
||||
PetscErrorCode :: ierr
|
||||
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
|
||||
|
|
@ -324,10 +326,10 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
|
|||
|
||||
Fdot = utilities_calculateRate(guess, &
|
||||
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
|
||||
math_rotate_backward33(F_aimDot,rotation_BC))
|
||||
rotation_BC%rotTensor2(F_aimDot,active=.true.))
|
||||
F_tauDot = utilities_calculateRate(guess, &
|
||||
F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
|
||||
math_rotate_backward33(F_aimDot,rotation_BC))
|
||||
rotation_BC%rotTensor2(F_aimDot,active=.true.))
|
||||
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
|
||||
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
|
||||
|
||||
|
|
@ -338,7 +340,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
|
|||
! update average and local deformation gradients
|
||||
F_aim = F_aim_lastInc + F_aimDot * timeinc
|
||||
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
|
||||
math_rotate_backward33(F_aim,rotation_BC)),&
|
||||
rotation_BC%rotTensor2(F_aim,active=.true.)),&
|
||||
[9,grid(1),grid(2),grid3])
|
||||
if (guess) then
|
||||
F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), &
|
||||
|
|
@ -349,8 +351,8 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
|
|||
F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, &
|
||||
math_mul3333xx33(C_scale,&
|
||||
matmul(transpose(F_lambda33),&
|
||||
F_lambda33)-math_I3))*0.5_pReal)&
|
||||
+ math_I3
|
||||
F_lambda33)-math_I3))*0.5_pReal) &
|
||||
+ math_I3
|
||||
F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
|
||||
enddo; enddo; enddo
|
||||
endif
|
||||
|
|
@ -514,7 +516,7 @@ subroutine formResidual(in, FandF_tau, &
|
|||
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
|
||||
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
||||
' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
|
||||
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
||||
' deformation gradient aim =', transpose(F_aim)
|
||||
flush(6)
|
||||
|
|
@ -533,7 +535,7 @@ subroutine formResidual(in, FandF_tau, &
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! doing convolution in Fourier space
|
||||
call utilities_FFTtensorForward
|
||||
call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC))
|
||||
call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(polarBeta*F_aim,active=.true.))
|
||||
call utilities_FFTtensorBackward
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -544,14 +546,14 @@ subroutine formResidual(in, FandF_tau, &
|
|||
! evaluate constitutive response
|
||||
call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory)
|
||||
P_av,C_volAvg,C_minMaxAvg, &
|
||||
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC)
|
||||
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC%asMatrix())
|
||||
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! stress BC handling
|
||||
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc
|
||||
err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim &
|
||||
-math_rotate_forward33(F_av,params%rotation_BC)) + &
|
||||
-params%rotation_BC%rotTensor2(F_av)) + &
|
||||
params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc
|
||||
! calculate divergence
|
||||
tensorField_real = 0.0_pReal
|
||||
|
|
|
|||
|
|
@ -104,7 +104,8 @@ module spectral_utilities
|
|||
end type tLoadCase
|
||||
|
||||
type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase
|
||||
real(pReal), dimension(3,3) :: stress_mask, stress_BC, rotation_BC
|
||||
real(pReal), dimension(3,3) :: stress_mask, stress_BC
|
||||
type(rotation) :: rotation_BC
|
||||
real(pReal) :: timeinc
|
||||
real(pReal) :: timeincOld
|
||||
end type tSolutionParams
|
||||
|
|
|
|||
14
src/math.f90
14
src/math.f90
|
|
@ -1315,20 +1315,6 @@ pure function math_rotate_forward33(tensor,R)
|
|||
end function math_rotate_forward33
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief rotate 33 tensor backward
|
||||
!> @details deprecated
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function math_rotate_backward33(tensor,R)
|
||||
|
||||
real(pReal), dimension(3,3) :: math_rotate_backward33
|
||||
real(pReal), dimension(3,3), intent(in) :: tensor, R
|
||||
|
||||
math_rotate_backward33 = matmul(transpose(R),matmul(tensor,R))
|
||||
|
||||
end function math_rotate_backward33
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief rotate 3333 tensor C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop
|
||||
!> @details deprecated
|
||||
|
|
|
|||
|
|
@ -163,7 +163,7 @@ subroutine fromQuaternion(self,qu)
|
|||
class(rotation), intent(out) :: self
|
||||
real(pReal), dimension(4), intent(in) :: qu
|
||||
|
||||
if (dNeq(norm2(qu),1.0)) &
|
||||
if (dNeq(norm2(qu),1.0_pReal)) &
|
||||
call IO_error(402,ext_msg='fromQuaternion')
|
||||
|
||||
self%q = qu
|
||||
|
|
|
|||
Loading…
Reference in New Issue