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

P_aim should be independent from P_av 

P_av is not defined after restart or cutback. Restart with change of
load case is probably still an issue
This commit is contained in:
Martin Diehl 2020-12-15 08:06:50 +01:00
parent b6db675f1a
commit 2dd520b4a2
8 changed files with 224 additions and 224 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 08f8aea465a1b5e476b584bcae7927d113919b1d
Subproject commit de65e1df5a76362de93667e9820dbf330b56f96d

4
python/damask/grid_filters.py
View File

@ -234,7 +234,7 @@ def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True):
origin[_np.where(cells==1)] = 0.0
if cells.prod() != len(coordinates0):
raise ValueError('Data count {len(coordinates0)} does not match cells {cells}.')
raise ValueError(f'Data count {len(coordinates0)} does not match cells {cells}.')
start = origin + delta*.5
end = origin - delta*.5 + size
@ -387,7 +387,7 @@ def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
origin = mincorner
if (cells+1).prod() != len(coordinates0):
raise ValueError('Data count {len(coordinates0)} does not match cells {cells}.')
raise ValueError(f'Data count {len(coordinates0)} does not match cells {cells}.')
atol = _np.max(size)*5e-2
if not (_np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],cells[0]+1),atol=atol) and \

111
src/grid/DAMASK_grid.f90
View File

@ -42,8 +42,8 @@ program DAMASK_grid
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors
real(pReal), dimension(9) :: temp_valueVector !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
logical, dimension(9) :: temp_maskVector !< temporarily from loadcase file when reading in tensors
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
@ -143,8 +143,6 @@ program DAMASK_grid
mech_restartWrite => grid_mech_spectral_basic_restartWrite
case ('Polarisation')
if(debug_grid%contains('basic')) &
call IO_warning(42, ext_msg='debug Divergence')
mech_init => grid_mech_spectral_polarisation_init
mech_forward => grid_mech_spectral_polarisation_forward
mech_solution => grid_mech_spectral_polarisation_solution
@ -152,8 +150,6 @@ program DAMASK_grid
mech_restartWrite => grid_mech_spectral_polarisation_restartWrite
case ('FEM')
if(debug_grid%contains('basic')) &
call IO_warning(42, ext_msg='debug Divergence')
mech_init => grid_mech_FEM_init
mech_forward => grid_mech_FEM_forward
mech_solution => grid_mech_FEM_solution
@ -178,11 +174,11 @@ program DAMASK_grid
allocate(loadCases(l)%ID(nActiveFields))
field = 1
loadCases(l)%ID(field) = FIELD_MECH_ID ! mechanical active by default
thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
field = field + 1
loadCases(l)%ID(field) = FIELD_THERMAL_ID
endif thermalActive
damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then
damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then
field = field + 1
loadCases(l)%ID(field) = FIELD_DAMAGE_ID
endif damageActive
@ -190,33 +186,35 @@ program DAMASK_grid
load_step => load_steps%get(l)
step_mech => load_step%get('mechanics')
loadCases(l)%stress%myType='P'
loadCases(l)%stress%myType=''
readMech: do m = 1, step_mech%length
select case (step_mech%getKey(m))
case('dot_F','L','F') ! assign values for the deformation BC matrix
case ('L','dot_F','F') ! assign values for the deformation BC matrix
loadCases(l)%deformation%myType = step_mech%getKey(m)
temp_valueVector = 0.0_pReal
step_deformation => step_mech%get(m)
do j = 1, 9
temp_maskVector(j) = step_deformation%get_asString(j) /= 'x' ! true if not a 'x'
if (temp_maskVector(j)) temp_valueVector(j) = step_deformation%get_asFloat(j) ! read value where applicable
enddo
loadCases(l)%deformation%mask = transpose(reshape(temp_maskVector,[ 3,3])) ! mask in 3x3 notation
loadCases(l)%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation
case('P')
temp_valueVector = 0.0_pReal
step_stress => step_mech%get(m)
do j = 1, 9
temp_maskVector(j) = step_stress%get_asString(j) /= 'x' ! true if not a 'x'
if (temp_maskVector(j)) temp_valueVector(j) = step_stress%get_asFloat(j) ! read value where applicable
temp_maskVector(j) = step_deformation%get_asString(j) /= 'x'
if (temp_maskVector(j)) temp_valueVector(j) = step_deformation%get_asFloat(j)
enddo
loadCases(l)%stress%mask = transpose(reshape(temp_maskVector,[ 3,3]))
loadCases(l)%deformation%mask = transpose(reshape(temp_maskVector,[3,3]))
loadCases(l)%deformation%values = math_9to33(temp_valueVector)
case ('dot_P','P')
loadCases(l)%stress%myType = step_mech%getKey(m)
step_stress => step_mech%get(m)
temp_valueVector = 0.0_pReal
do j = 1, 9
temp_maskVector(j) = step_stress%get_asString(j) /= 'x'
if (temp_maskVector(j)) temp_valueVector(j) = step_stress%get_asFloat(j)
enddo
loadCases(l)%stress%mask = transpose(reshape(temp_maskVector,[3,3]))
loadCases(l)%stress%values = math_9to33(temp_valueVector)
end select
call loadCases(l)%rot%fromAxisAngle(step_mech%get_asFloats('R',defaultVal = real([0.0,0.0,1.0,0.0],pReal)),degrees=.true.)
enddo readMech
if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/F/dot_F missing')
if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/dot_F/F missing')
step_discretization => load_step%get('discretization')
if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing')
@ -239,50 +237,60 @@ program DAMASK_grid
if (any(loadCases(l)%deformation%mask(j,1:3) .eqv. .true.) .and. &
any(loadCases(l)%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined
enddo
print*, ' L:'
else if (loadCases(l)%deformation%myType == 'F') then
endif
if (loadCases(l)%deformation%myType == 'F') then
print*, ' F:'
else if (loadCases(l)%deformation%myType == 'dot_F') then
print*, ' dot_F:'
else
print*, ' '//loadCases(l)%deformation%myType//' / 1/s:'
endif
do i = 1, 3; do j = 1, 3
if(loadCases(l)%deformation%mask(i,j)) then
if (loadCases(l)%deformation%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j)
else
write(IO_STDOUT,'(2x,12a)',advance='no') ' x '
endif
enddo; write(IO_STDOUT,'(/)',advance='no')
enddo
if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831 ! exclusive or masking only
if (any(loadCases(l)%stress%mask .and. transpose(loadCases(l)%stress%mask) .and. (math_I3<1))) &
errorID = 838 ! no rotation is allowed by stress BC
print*, ' P / MPa:'
do i = 1, 3; do j = 1, 3
if(loadCases(l)%stress%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal
else
write(IO_STDOUT,'(2x,12a)',advance='no') ' x '
endif
enddo; write(IO_STDOUT,'(/)',advance='no')
enddo
if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831
if (any(loadCases(l)%stress%mask .and. transpose(loadCases(l)%stress%mask) .and. (math_I3<1))) &
errorID = 838 ! no rotation is allowed by stress BC
if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:'
if (loadCases(l)%stress%myType == 'dot_P') print*, ' dot_P / MPa/s:'
if (loadCases(l)%stress%myType /= '') then
do i = 1, 3; do j = 1, 3
if (loadCases(l)%stress%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal
else
write(IO_STDOUT,'(2x,12a)',advance='no') ' x '
endif
enddo; write(IO_STDOUT,'(/)',advance='no')
enddo
endif
if (any(dNeq(loadCases(l)%rot%asMatrix(), math_I3))) &
write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'R:',&
transpose(loadCases(l)%rot%asMatrix())
if (loadCases(l)%t < 0.0_pReal) errorID = 834
print'(a,f0.3)', ' t: ', loadCases(l)%t
if (loadCases(l)%N < 1) errorID = 835
print'(a,i0)', ' N: ', loadCases(l)%N
if (loadCases(l)%f_out < 1) errorID = 836
print'(a,i0)', ' f_out: ', loadCases(l)%f_out
if (loadCases(l)%r <= 0.0) errorID = 833
print'(a,f0.3)', ' r: ', loadCases(l)%r
if (loadCases(l)%r <= 0.0) errorID = 833
if (loadCases(l)%t < 0.0_pReal) errorID = 834
if (loadCases(l)%N < 1) errorID = 835
if (loadCases(l)%f_out < 1) errorID = 836
if (loadCases(l)%f_restart < 1) errorID = 839
if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then
print'(a)', ' r: 1 (constant step widths)'
else
print'(a,f0.3)', ' r: ', loadCases(l)%r
endif
print'(a,f0.3)', ' t: ', loadCases(l)%t
print'(a,i0)', ' N: ', loadCases(l)%N
print'(a,i0)', ' f_out: ', loadCases(l)%f_out
if (loadCases(l)%f_restart < huge(0)) &
print'(a,i0)', ' f_restart: ', loadCases(l)%f_restart
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
endif reportAndCheck
enddo
@ -309,8 +317,6 @@ program DAMASK_grid
writeHeader: if (interface_restartInc < 1) then
open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
if (debug_grid%contains('basic')) print'(/,a)', ' header of statistics file written out'
flush(IO_STDOUT)
else writeHeader
open(newunit=statUnit,file=trim(getSolverJobName())//&
'.sta',form='FORMATTED', position='APPEND', status='OLD')
@ -319,6 +325,7 @@ program DAMASK_grid
writeUndeformed: if (interface_restartInc < 1) then
print'(/,a)', ' ... writing initial configuration to file ........................'
flush(IO_STDOUT)
call CPFEM_results(0,0.0_pReal)
endif writeUndeformed

123
src/grid/grid_mech_FEM.f90
View File

@ -31,16 +31,16 @@ module grid_mech_FEM
type :: tNumerics
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
real(pReal) :: &
eps_div_atol, & !< absolute tolerance for equilibrium
eps_div_rtol, & !< relative tolerance for equilibrium
eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC
eps_stress_rtol !< relative tolerance for fullfillment of stress BC
eps_div_atol, & !< absolute tolerance for equilibrium
eps_div_rtol, & !< relative tolerance for equilibrium
eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC
eps_stress_rtol !< relative tolerance for fullfillment of stress BC
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
type(tNumerics) :: num ! numerics parameters. Better name?
logical :: debugRotation
@ -64,7 +64,7 @@ module grid_mech_FEM
real(pReal), dimension(3,3) :: &
F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
P_aim = 0.0_pReal
character(len=:), allocatable :: incInfo !< time and increment information
@ -93,10 +93,8 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_init
real(pReal) :: HGCoeff = 0.0e-2_pReal
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
real(pReal), dimension(4,8) :: &
real(pReal), parameter :: HGCoeff = 0.0e-2_pReal
real(pReal), parameter, dimension(4,8) :: &
HGcomp = reshape([ 1.0_pReal, 1.0_pReal, 1.0_pReal,-1.0_pReal, &
1.0_pReal,-1.0_pReal,-1.0_pReal, 1.0_pReal, &
-1.0_pReal, 1.0_pReal,-1.0_pReal, 1.0_pReal, &
@ -121,18 +119,19 @@ subroutine grid_mech_FEM_init
!-------------------------------------------------------------------------------------------------
! debugging options
debug_grid => config_debug%get('grid', defaultVal=emptyList)
debug_grid => config_debug%get('grid',defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
!-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal)
num%itmin = num_grid%get_asInt ('itmin', defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
@ -225,6 +224,7 @@ subroutine grid_mech_FEM_init
fileHandle = HDF5_openFile(fileName)
groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_read(groupHandle,P_aim, 'P_aim')
call HDF5_read(groupHandle,F_aim, 'F_aim')
call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(groupHandle,F_aimDot, 'F_aimDot')
@ -238,9 +238,9 @@ subroutine grid_mech_FEM_init
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(F)
call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F
0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr)
@ -295,6 +295,7 @@ function grid_mech_FEM_solution(incInfoIn) result(solution)
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask)
end function grid_mech_FEM_solution
@ -302,34 +303,26 @@ end function grid_mech_FEM_solution
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: &
logical, intent(in) :: &
cutBack, &
guess
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
real(pReal), intent(in) :: &
Delta_t_old, &
Delta_t, &
t_remaining !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(rotation), intent(in) :: &
type(rotation), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
@ -339,7 +332,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
else
C_volAvgLastInc = C_volAvg
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components
F_aim_lastInc = F_aim
!-----------------------------------------------------------------------------------------------
@ -347,18 +340,18 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F
F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
elseif(deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot &
elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
+ merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask)
endif
if (guess) then
call VecWAXPY(solution_rate,-1.0_pReal,solution_lastInc,solution_current,ierr)
CHKERRQ(ierr)
call VecScale(solution_rate,1.0_pReal/timeinc_old,ierr); CHKERRQ(ierr)
call VecScale(solution_rate,1.0_pReal/Delta_t_old,ierr); CHKERRQ(ierr)
else
call VecSet(solution_rate,0.0_pReal,ierr); CHKERRQ(ierr)
endif
@ -371,23 +364,28 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
if (stress_BC%myType=='P') then
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
elseif (stress_BC%myType=='dot_P') then !UNTESTED
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
endif
F_aim = F_aim_lastInc + F_aimDot * Delta_t
if (stress_BC%myType=='P') P_aim = P_aim &
+ merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t
if (stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t
call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr)
call VecAXPY(solution_current,Delta_t,solution_rate,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = Delta_t
end subroutine grid_mech_FEM_forward
!--------------------------------------------------------------------------------------------------
!> @brief Age
!> @brief Update coordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_updateCoords
@ -415,6 +413,7 @@ subroutine grid_mech_FEM_restartWrite
fileHandle = HDF5_openFile(fileName,'w')
groupHandle = HDF5_addGroup(fileHandle,'solver')
call HDF5_write(groupHandle,P_aim, 'P_aim')
call HDF5_write(groupHandle,F_aim, 'F_aim')
call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_write(groupHandle,F_aimDot, 'F_aimDot')
@ -441,11 +440,11 @@ end subroutine grid_mech_FEM_restartWrite
subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,ierr)
SNES :: snes_local
PetscInt, intent(in) :: PETScIter
PetscReal, intent(in) :: &
devNull1, &
devNull2, &
fnorm
PetscInt, intent(in) :: PETScIter
PetscReal, intent(in) :: &
devNull1, &
devNull2, &
fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
@ -458,10 +457,10 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol)
BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
if ((totalIter >= num%itmin .and. &
all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
if (terminallyIll .or. &
(totalIter >= num%itmin .and. &
all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal))) then
reason = 1
elseif (totalIter >= num%itmax) then
reason = -1
@ -666,16 +665,16 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ
call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr)
CHKERRQ(ierr)
call DMGetGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr)
call DMGetGlobalVector(da_local,coordinates,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr); CHKERRQ(ierr)
ele = 0
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
ele = ele + 1
x_scal(0:2,i,j,k) = discretization_IPcoords(1:3,ele)
enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates)
call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes
call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr); CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates)
call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr); CHKERRQ(ierr) ! get rigid body deformation modes
call DMRestoreGlobalVector(da_local,coordinates,ierr); CHKERRQ(ierr)
call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)

97
src/grid/grid_mech_spectral_basic.f90
View File

@ -94,8 +94,6 @@ contains
subroutine grid_mech_spectral_basic_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
F ! pointer to solution data
@ -118,20 +116,20 @@ subroutine grid_mech_spectral_basic_init
!-------------------------------------------------------------------------------------------------
! debugging options
debug_grid => config_debug%get('grid', defaultVal=emptyList)
debug_grid => config_debug%get('grid',defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
!-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=1.0e-3_pReal)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
@ -139,7 +137,7 @@ subroutine grid_mech_spectral_basic_init
if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol')
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin')
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
@ -149,14 +147,14 @@ subroutine grid_mech_spectral_basic_init
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
localK = 0
localK = 0
localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
@ -189,6 +187,7 @@ subroutine grid_mech_spectral_basic_init
fileHandle = HDF5_openFile(fileName)
groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_read(groupHandle,P_aim, 'P_aim')
call HDF5_read(groupHandle,F_aim, 'F_aim')
call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(groupHandle,F_aimDot, 'F_aimDot')
@ -200,9 +199,9 @@ subroutine grid_mech_spectral_basic_init
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer
@ -262,6 +261,7 @@ function grid_mech_spectral_basic_solution(incInfoIn) result(solution)
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask)
end function grid_mech_spectral_basic_solution
@ -269,32 +269,25 @@ end function grid_mech_spectral_basic_solution
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: &
logical, intent(in) :: &
cutBack, &
guess
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
real(pReal), intent(in) :: &
Delta_t_old, &
Delta_t, &
t_remaining !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(rotation), intent(in) :: &
type(rotation), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: F
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
@ -305,7 +298,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components
F_aim_lastInc = F_aim
!-----------------------------------------------------------------------------------------------
@ -313,40 +306,45 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F
F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
elseif(deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot &
elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
+ merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask)
endif
Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
rotation_BC%rotate(F_aimDot,active=.true.))
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3])
endif
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
if (stress_BC%myType=='P') then
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
elseif (stress_BC%myType=='dot_P') then !UNTESTED
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
endif
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
F_aim = F_aim_lastInc + F_aimDot * Delta_t
if (stress_BC%myType=='P') P_aim = P_aim &
+ merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t
if (stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = Delta_t
end subroutine grid_mech_spectral_basic_forward
!--------------------------------------------------------------------------------------------------
!> @brief Age
!> @brief Update coordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_updateCoords
@ -378,6 +376,7 @@ subroutine grid_mech_spectral_basic_restartWrite
fileHandle = HDF5_openFile(fileName,'w')
groupHandle = HDF5_addGroup(fileHandle,'solver')
call HDF5_write(groupHandle,P_aim, 'P_aim')
call HDF5_write(groupHandle,F_aim, 'F_aim')
call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_write(groupHandle,F_aimDot, 'F_aimDot')
@ -463,7 +462,7 @@ subroutine formResidual(in, F, &
PetscErrorCode :: ierr
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
@ -472,7 +471,7 @@ subroutine formResidual(in, F, &
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if (debugRotation) &
if(debugRotation) &
write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &

107
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -105,8 +105,6 @@ contains
subroutine grid_mech_spectral_polarisation_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
FandF_tau, & ! overall pointer to solution data
@ -147,16 +145,16 @@ subroutine grid_mech_spectral_polarisation_init
num%alpha = num_grid%get_asFloat('alpha', defaultVal=1.0_pReal)
num%beta = num_grid%get_asFloat('beta', defaultVal=1.0_pReal)
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
if (num%eps_curl_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_curl_atol')
if (num%eps_curl_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_curl_rtol')
if (num%eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol')
if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol')
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin')
if (num%alpha <= 0.0_pReal .or. num%alpha > 2.0_pReal) call IO_error(301,ext_msg='alpha')
if (num%beta < 0.0_pReal .or. num%beta > 2.0_pReal) call IO_error(301,ext_msg='beta')
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
if (num%eps_curl_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_curl_atol')
if (num%eps_curl_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_curl_rtol')
if (num%eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol')
if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol')
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin')
if (num%alpha <= 0.0_pReal .or. num%alpha > 2.0_pReal) call IO_error(301,ext_msg='alpha')
if (num%beta < 0.0_pReal .or. num%beta > 2.0_pReal) call IO_error(301,ext_msg='beta')
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
@ -211,6 +209,7 @@ subroutine grid_mech_spectral_polarisation_init
fileHandle = HDF5_openFile(fileName)
groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_read(groupHandle,P_aim, 'P_aim')
call HDF5_read(groupHandle,F_aim, 'F_aim')
call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(groupHandle,F_aimDot, 'F_aimDot')
@ -226,9 +225,9 @@ subroutine grid_mech_spectral_polarisation_init
F_tau_lastInc = 2.0_pReal*F_lastInc
endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
@ -294,6 +293,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution)
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask)
end function grid_mech_spectral_polarisation_solution
@ -301,34 +301,27 @@ end function grid_mech_spectral_polarisation_solution
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: &
logical, intent(in) :: &
cutBack, &
guess
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
real(pReal), intent(in) :: &
Delta_t_old, &
Delta_t, &
t_remaining !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
type(rotation), intent(in) :: &
type(rotation), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau
integer :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
F => FandF_tau(0: 8,:,:,:)
@ -341,7 +334,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components
F_aim_lastInc = F_aim
!-----------------------------------------------------------------------------------------------
@ -349,19 +342,19 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F
F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
elseif(deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot &
elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
+ merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask)
endif
Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
rotation_BC%rotate(F_aimDot,active=.true.))
F_tauDot = utilities_calculateRate(guess, &
F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), Delta_t_old, &
rotation_BC%rotate(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])
@ -371,38 +364,41 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
if (stress_BC%myType=='P') then
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
elseif (stress_BC%myType=='dot_P') then !UNTESTED
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
endif
F_aim = F_aim_lastInc + F_aimDot * Delta_t
if(stress_BC%myType=='P') P_aim = P_aim &
+ merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t
if(stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),&
[9,grid(1),grid(2),grid3])
if (guess) then
F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), &
F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), &
[9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition
else
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
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 &
+ math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, &
math_mul3333xx33(C_scale,matmul(transpose(F_lambda33),F_lambda33)-math_I3)))
F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
enddo; enddo; enddo
endif
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = Delta_t
end subroutine grid_mech_spectral_polarisation_forward
!--------------------------------------------------------------------------------------------------
!> @brief Age
!> @brief Update coordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_updateCoords
@ -436,6 +432,7 @@ subroutine grid_mech_spectral_polarisation_restartWrite
fileHandle = HDF5_openFile(fileName,'w')
groupHandle = HDF5_addGroup(fileHandle,'solver')
call HDF5_write(groupHandle,F_aim, 'P_aim')
call HDF5_write(groupHandle,F_aim, 'F_aim')
call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_write(groupHandle,F_aimDot, 'F_aimDot')
@ -480,11 +477,11 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
divTol = max(maxval(abs(P_av)) *num%eps_div_rtol ,num%eps_div_atol)
BCTol = max(maxval(abs(P_av)) *num%eps_stress_rtol,num%eps_stress_atol)
if ((totalIter >= num%itmin .and. &
all([ err_div /divTol, &
err_curl/curlTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
if (terminallyIll .or. &
(totalIter >= num%itmin .and. &
all([ err_div /divTol, &
err_curl/curlTol, &
err_BC /BCTol ] < 1.0_pReal))) then
reason = 1
elseif (totalIter >= num%itmax) then
reason = -1
@ -555,7 +552,7 @@ subroutine formResidual(in, FandF_tau, &
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if(debugRotation) &
if (debugRotation) &
write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &

2
src/grid/spectral_utilities.f90
View File

@ -93,7 +93,7 @@ module spectral_utilities
real(pReal), dimension(3,3) :: stress_BC
logical, dimension(3,3) :: stress_mask
type(rotation) :: rotation_BC
real(pReal) :: timeinc, timeincOld
real(pReal) :: timeinc
end type tSolutionParams
type :: tNumerics

2
src/mesh/mesh_mech_FEM.f90
View File

@ -32,7 +32,6 @@ module mesh_mech_FEM
type tSolutionParams
type(tFieldBC) :: fieldBC
real(pReal) :: timeinc
real(pReal) :: timeincOld
end type tSolutionParams
type(tSolutionParams) :: params
@ -302,7 +301,6 @@ type(tSolutionState) function FEM_mech_solution( &
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%timeinc = timeinc
params%timeincOld = timeinc_old
params%fieldBC = fieldBC
call SNESSolve(mech_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mech_snes based on solution guess (result in solution)