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

num structure replicated; polishing

This commit is contained in:
Sharan Roongta 2020-06-24 16:45:13 +02:00
parent 434bfffc46
commit be84561e2e
9 changed files with 190 additions and 239 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/crystallite.f90
View File

@ -83,7 +83,7 @@ module crystallite
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
nState, & !< state loop limit nState, & !< state loop limit
nStress !< stress loop limit nStress !< stress loop limit
character(len=pStringLen) :: & character(len=:), allocatable :: &
integrator !< integrator scheme integrator !< integrator scheme
real(pReal) :: & real(pReal) :: &
subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback
@ -196,7 +196,7 @@ subroutine crystallite_init
if(num%nStress< 1) call IO_error(301,ext_msg='nStress') if(num%nStress< 1) call IO_error(301,ext_msg='nStress')
select case(trim(num%integrator)) select case(num%integrator)
case('FPI') case('FPI')
integrateState => integrateStateFPI integrateState => integrateStateFPI
case('Euler') case('Euler')

91
src/grid/grid_mech_FEM.f90
View File

@ -30,6 +30,23 @@ module grid_mech_FEM
! derived types ! derived types
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
type, private :: tNumerics
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
real(pReal) :: &
err_div, &
divTol, &
BCTol, &
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
character(len=:), allocatable :: &
petsc_options
end type tNumerics
type(tNumerics), private :: num
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM, private :: mech_grid DM, private :: mech_grid
@ -97,8 +114,7 @@ subroutine grid_mech_FEM_init
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: & character(len=pStringLen) :: &
fileName, & fileName
petsc_options
class(tNode), pointer :: & class(tNode), pointer :: &
num_grid num_grid
real(pReal), dimension(3,3,3,3) :: devNull real(pReal), dimension(3,3,3,3) :: devNull
@ -108,16 +124,30 @@ subroutine grid_mech_FEM_init
write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>'; flush(6)
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameter ! read numerical parameter and do sanity checks
num_grid => numerics_root%get('grid',defaultVal=emptyDict) num_grid => numerics_root%get('grid',defaultVal=emptyDict)
petsc_options = num_grid%get_asString('petsc_options',defaultVal='') num%petsc_options = num_grid%get_asString('petsc_options', defaultVal='')
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=0.01_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')
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')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres &
&-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) &-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -419,49 +449,23 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
SNESConvergedReason :: reason SNESConvergedReason :: reason
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
real(pReal) :: & real(pReal) :: &
err_div, & err_div, &
divTol, & divTol, &
BCTol, & BCTol
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
class(tNode), pointer :: &
num_grid
!-----------------------------------------------------------------------------------
! reading numerical parameters and do sanity check
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
eps_div_atol = num_grid%get_asFloat('eps_div_atol',defaultVal=1.0e-4_pReal)
eps_div_rtol = num_grid%get_asFloat('eps_div_rtol',defaultVal=5.0e-4_pReal)
eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal)
itmin = num_grid%get_asInt('itmin',defaultVal=1)
itmax = num_grid%get_asInt('itmax',defaultVal=250)
if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
if (eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol')
if (eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol')
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
!------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------
err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ
divTol = max(maxval(abs(P_av))*eps_div_rtol ,eps_div_atol) divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol)
BCTol = max(maxval(abs(P_av))*eps_stress_rtol,eps_stress_atol) BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
if ((totalIter >= itmin .and. & if ((totalIter >= num%itmin .and. &
all([ err_div/divTol, & all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal)) & err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then .or. terminallyIll) then
reason = 1 reason = 1
elseif (totalIter >= itmax) then elseif (totalIter >= num%itmax) then
reason = -1 reason = -1
else else
reason = 0 reason = 0
@ -499,19 +503,6 @@ subroutine formResidual(da_local,x_local, &
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
real(pReal), dimension(3,3,3,3) :: devNull real(pReal), dimension(3,3,3,3) :: devNull
integer :: &
itmin, &
itmax
class(tNode), pointer :: &
num_grid
!----------------------------------------------------------------------
! read numerical paramteters and do sanity checks
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
itmin = num_grid%get_asInt('itmin',defaultVal=1)
itmax = num_grid%get_asInt('itmax',defaultVal=250)
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr)
@ -522,7 +513,7 @@ subroutine formResidual(da_local,x_local, &
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))

91
src/grid/grid_mech_spectral_basic.f90
View File

@ -32,9 +32,19 @@ module grid_mech_spectral_basic
type, private :: tNumerics type, private :: tNumerics
logical :: update_gamma !< update gamma operator with current stiffness logical :: update_gamma !< update gamma operator with current stiffness
integer :: &
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
character(len=:), allocatable :: &
petsc_options
end type tNumerics end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics), private :: num ! numerics parameters. Better name?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
@ -98,8 +108,7 @@ subroutine grid_mech_spectral_basic_init
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
integer :: fileUnit integer :: fileUnit
character(len=pStringLen) :: & character(len=pStringLen) :: &
fileName, & fileName
petsc_options
write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(6)
@ -110,16 +119,31 @@ subroutine grid_mech_spectral_basic_init
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters ! read numerical parameters and do sanity checks
num_grid => numerics_root%get('grid',defaultVal=emptyDict) num_grid => numerics_root%get('grid',defaultVal=emptyDict)
petsc_options = num_grid%get_asString('petsc_options',defaultVal='')
num%update_gamma = num_grid%get_asBool('update_gamma',defaultVal=.false.) num%petsc_options = num_grid%get_asString ('petsc_options', defaultVal='')
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=0.01_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')
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')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -389,47 +413,19 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
SNESConvergedReason :: reason SNESConvergedReason :: reason
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
real(pReal) :: & real(pReal) :: &
divTol, & divTol, &
BCTol, & BCTol
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
class(tNode), pointer :: &
num_grid
!----------------------------------------------------------------------------------- divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol)
! reading numerical parameters and do sanity check BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
eps_div_atol = num_grid%get_asFloat('eps_div_atol',defaultVal=1.0e-4_pReal)
eps_div_rtol = num_grid%get_asFloat('eps_div_rtol',defaultVal=5.0e-4_pReal)
eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal)
itmin = num_grid%get_asInt('itmin',defaultVal=1) if ((totalIter >= num%itmin .and. &
itmax = num_grid%get_asInt('itmax',defaultVal=250)
if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
if (eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol')
if (eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol')
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
!------------------------------------------------------------------------------------
divTol = max(maxval(abs(P_av))*eps_div_rtol ,eps_div_atol)
BCTol = max(maxval(abs(P_av))*eps_stress_rtol,eps_stress_atol)
if ((totalIter >= itmin .and. &
all([ err_div/divTol, & all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal)) & err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then .or. terminallyIll) then
reason = 1 reason = 1
elseif (totalIter >= itmax) then elseif (totalIter >= num%itmax) then
reason = -1 reason = -1
else else
reason = 0 reason = 0
@ -466,19 +462,6 @@ subroutine formResidual(in, F, &
nfuncs nfuncs
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: &
itmin, &
itmax
class(tNode), pointer :: &
num_grid
!----------------------------------------------------------------------
! read numerical paramteter and do sanity checks
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
itmin = num_grid%get_asInt('itmin',defaultVal=1)
itmax = num_grid%get_asInt('itmax',defaultVal=250)
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
@ -488,7 +471,7 @@ subroutine formResidual(in, F, &
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))

128
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -33,9 +33,24 @@ module grid_mech_spectral_polarisation
type, private :: tNumerics type, private :: tNumerics
logical :: update_gamma !< update gamma operator with current stiffness logical :: update_gamma !< update gamma operator with current stiffness
character(len=:), allocatable :: &
petsc_options
integer :: &
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_curl_atol, & !< absolute tolerance for compatibility
eps_curl_rtol, & !< relative tolerance for compatibility
eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC
eps_stress_rtol !< relative tolerance for fullfillment of stress BC
real(pReal) :: &
polarAlpha, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
polarBeta !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
end type tNumerics end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics), private :: num ! numerics parameters. Better name?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
@ -107,8 +122,7 @@ subroutine grid_mech_spectral_polarisation_init
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
integer :: fileUnit integer :: fileUnit
character(len=pStringLen) :: & character(len=pStringLen) :: &
fileName, & fileName
petsc_options
write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(6)
@ -118,14 +132,36 @@ subroutine grid_mech_spectral_polarisation_init
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters ! read numerical parameters
num_grid => numerics_root%get('grid',defaultVal=emptyDict) num_grid => numerics_root%get('grid',defaultVal=emptyDict)
petsc_options = num_grid%get_asString('petsc_options',defaultVal='')
num%update_gamma = num_grid%get_asBool('update_gamma',defaultVal=.false.) num%petsc_options = num_grid%get_asString('petsc_options', defaultVal='')
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_curl_atol = num_grid%get_asFloat ('eps_curl_atol', defaultVal=1.0e-10_pReal)
num%eps_curl_rtol = num_grid%get_asFloat ('eps_curl_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=0.01_pReal)
num%itmin = num_grid%get_asInt ('itmin', defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
num%polarAlpha = num_grid%get_asFloat ('polaralpha', defaultVal=1.0_pReal)
num%polarBeta = num_grid%get_asFloat ('polarbeta', 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%polarAlpha <= 0.0_pReal .or. num%polarAlpha > 2.0_pReal) call IO_error(301,ext_msg='polarAlpha')
if (num%polarBeta < 0.0_pReal .or. num%polarBeta > 2.0_pReal) call IO_error(301,ext_msg='polarBeta')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -436,56 +472,22 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
SNESConvergedReason :: reason SNESConvergedReason :: reason
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
real(pReal) :: & real(pReal) :: &
curlTol, & curlTol, &
divTol, & divTol, &
BCTol, & BCTol
eps_div_atol, & !< absolute tolerance for equilibrium
eps_div_rtol, & !< relative tolerance for equilibrium
eps_curl_atol, & !< absolute tolerance for compatibility
eps_curl_rtol, & !< relative tolerance for compatibility
eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC
eps_stress_rtol !< relative tolerance for fullfillment of stress BC
class(tNode), pointer :: &
num_grid
!----------------------------------------------------------------------------------- curlTol = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol ,num%eps_curl_atol)
! reading numerical parameters and do sanity check divTol = max(maxval(abs(P_av)) *num%eps_div_rtol ,num%eps_div_atol)
num_grid => numerics_root%get('grid',defaultVal=emptyDict) BCTol = max(maxval(abs(P_av)) *num%eps_stress_rtol,num%eps_stress_atol)
eps_div_atol = num_grid%get_asFloat('eps_div_atol',defaultVal=1.0e-4_pReal)
eps_div_rtol = num_grid%get_asFloat('eps_div_rtol',defaultVal=5.0e-4_pReal)
eps_curl_atol = num_grid%get_asFloat('eps_curl_atol',defaultVal=1.0e-10_pReal)
eps_curl_rtol = num_grid%get_asFloat('eps_curl_rtol',defaultVal=5.0e-4_pReal)
eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal)
itmin = num_grid%get_asInt('itmin',defaultVal=1) if ((totalIter >= num%itmin .and. &
itmax = num_grid%get_asInt('itmax',defaultVal=250)
if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
if (eps_curl_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_curl_atol')
if (eps_curl_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_curl_rtol')
if (eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol')
if (eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol')
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
!------------------------------------------------------------------------------------
curlTol = max(maxval(abs(F_aim-math_I3))*eps_curl_rtol ,eps_curl_atol)
divTol = max(maxval(abs(P_av)) *eps_div_rtol ,eps_div_atol)
BCTol = max(maxval(abs(P_av)) *eps_stress_rtol,eps_stress_atol)
if ((totalIter >= itmin .and. &
all([ err_div /divTol, & all([ err_div /divTol, &
err_curl/curlTol, & err_curl/curlTol, &
err_BC /BCTol ] < 1.0_pReal)) & err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then .or. terminallyIll) then
reason = 1 reason = 1
elseif (totalIter >= itmax) then elseif (totalIter >= num%itmax) then
reason = -1 reason = -1
else else
reason = 0 reason = 0
@ -527,27 +529,9 @@ subroutine formResidual(in, FandF_tau, &
nfuncs nfuncs
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
class(tNode), pointer :: & integer :: &
num_grid i, j, k, e
real(pReal) :: &
polarAlpha, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
polarBeta !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
integer :: &
i, j, k, e, &
itmin, itmax
!--------------------------------------------------------------------------------------------------
! read numerical paramteters and do sanity checks
num_grid => numerics_root%get('grid',defaultVal = emptyDict)
polarAlpha = num_grid%get_asFloat('polaralpha',defaultVal=1.0_pReal)
polarBeta = num_grid%get_asFloat('polarbeta', defaultVal=1.0_pReal)
itmin = num_grid%get_asInt('itmin',defaultVal=1)
itmax = num_grid%get_asInt('itmax',defaultVal=250)
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
if (polarAlpha <= 0.0_pReal .or. polarAlpha > 2.0_pReal) call IO_error(301,ext_msg='polarAlpha')
if (polarBeta < 0.0_pReal .or. polarBeta > 2.0_pReal) call IO_error(301,ext_msg='polarBeta')
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
F => FandF_tau(1:3,1:3,1,& F => FandF_tau(1:3,1:3,1,&
@ -570,7 +554,7 @@ subroutine formResidual(in, FandF_tau, &
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
@ -584,26 +568,26 @@ subroutine formResidual(in, FandF_tau, &
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
tensorField_real(1:3,1:3,i,j,k) = & tensorField_real(1:3,1:3,i,j,k) = &
polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& num%polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
polarAlpha*matmul(F(1:3,1:3,i,j,k), & num%polarAlpha*matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
enddo; enddo; enddo enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space ! doing convolution in Fourier space
call utilities_FFTtensorForward call utilities_FFTtensorForward
call utilities_fourierGammaConvolution(params%rotation_BC%rotate(polarBeta*F_aim,active=.true.)) call utilities_fourierGammaConvolution(params%rotation_BC%rotate(num%polarBeta*F_aim,active=.true.))
call utilities_FFTtensorBackward call utilities_FFTtensorBackward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) residual_F_tau = num%polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, & P_av,C_volAvg,C_minMaxAvg, &
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) F - residual_F_tau/num%polarBeta,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

10
src/grid/spectral_utilities.f90
View File

@ -121,10 +121,10 @@ module spectral_utilities
character(len=:), allocatable :: & character(len=:), allocatable :: &
spectral_derivative, & !< approximation used for derivatives in Fourier space spectral_derivative, & !< approximation used for derivatives in Fourier space
FFTW_plan_mode, & !< FFTW plan mode, see www.fftw.org FFTW_plan_mode, & !< FFTW plan mode, see www.fftw.org
PETSc_options petsc_options
end type tNumerics end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics), private :: num ! numerics parameters. Better name?
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
DERIVATIVE_CONTINUOUS_ID, & DERIVATIVE_CONTINUOUS_ID, &
@ -189,8 +189,6 @@ subroutine spectral_utilities_init
scalarSize = 1_C_INTPTR_T, & scalarSize = 1_C_INTPTR_T, &
vecSize = 3_C_INTPTR_T, & vecSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T tensorSize = 9_C_INTPTR_T
character(len=pStringLen) :: &
petsc_options
class(tNode), pointer :: & class(tNode), pointer :: &
num_grid num_grid
@ -220,13 +218,13 @@ subroutine spectral_utilities_init
' add more using the PETSc_Options keyword in numerics.config '; flush(6) ' add more using the PETSc_Options keyword in numerics.config '; flush(6)
num_grid => numerics_root%get('grid',defaultVal=emptyDict) num_grid => numerics_root%get('grid',defaultVal=emptyDict)
petsc_options = num_grid%get_asString('petsc_options',defaultVal='') num%petsc_options = num_grid%get_asString('petsc_options',defaultVal='')
call PETScOptionsClear(PETSC_NULL_OPTIONS,ierr) call PETScOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
if(debugPETSc) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) if(debugPETSc) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
grid1Red = grid(1)/2 + 1 grid1Red = grid(1)/2 + 1

4
src/marc/discretization_marc.f90
View File

@ -34,8 +34,8 @@ module discretization_marc
mesh_unitlength !< physical length of one unit in mesh mesh_unitlength !< physical length of one unit in mesh
integer, dimension(:), allocatable, public :: & integer, dimension(:), allocatable, public :: &
mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID
mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID
public :: & public :: &
discretization_marc_init discretization_marc_init

22
src/mesh/DAMASK_mesh.f90
View File

@ -49,7 +49,6 @@ program DAMASK_mesh
i, & i, &
errorID, & errorID, &
cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
maxCutBack, & !< max number of cutbacks
stepFraction = 0, & !< fraction of current time interval stepFraction = 0, & !< fraction of current time interval
currentLoadcase = 0, & !< current load case currentLoadcase = 0, & !< current load case
currentFace = 0, & currentFace = 0, &
@ -57,7 +56,6 @@ program DAMASK_mesh
totalIncsCounter = 0, & !< total # of increments totalIncsCounter = 0, & !< total # of increments
statUnit = 0, & !< file unit for statistics output statUnit = 0, & !< file unit for statistics output
stagIter, & stagIter, &
stagItMax, & !< max number of field level staggered iterations
component component
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
@ -65,6 +63,14 @@ program DAMASK_mesh
character(len=pStringLen) :: & character(len=pStringLen) :: &
incInfo, & incInfo, &
loadcase_string loadcase_string
type :: tNumerics
integer :: &
stagItMax, & !< max number of field level staggered iterations
maxCutBack !< max number of cutbacks
end type tNumerics
type(tNumerics) :: num
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet, field, dimPlex PetscInt :: faceSet, currentFaceSet, field, dimPlex
@ -81,11 +87,11 @@ program DAMASK_mesh
!--------------------------------------------------------------------- !---------------------------------------------------------------------
! reading field information from numerics file and do sanity checks ! reading field information from numerics file and do sanity checks
num_mesh => numerics_root%get('mesh', defaultVal=emptyDict) num_mesh => numerics_root%get('mesh', defaultVal=emptyDict)
stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10) num%stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10)
maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3) num%maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3)
if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (num%stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter')
if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') if (num%maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack')
! reading basic information from load case file and allocate data structure containing load cases ! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRA(ierr) !< dimension of mesh (2D or 3D) call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRA(ierr) !< dimension of mesh (2D or 3D)
@ -327,7 +333,7 @@ program DAMASK_mesh
enddo enddo
stagIter = stagIter + 1 stagIter = stagIter + 1
stagIterate = stagIter < stagItMax & stagIterate = stagIter < num%stagItMax &
.and. all(solres(:)%converged) & .and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
enddo enddo
@ -335,7 +341,7 @@ program DAMASK_mesh
! check solution ! check solution
cutBack = .False. cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back if (cutBackLevel < num%maxCutBack) then ! do cut back
write(6,'(/,a)') ' cut back detected' write(6,'(/,a)') ' cut back detected'
cutBack = .True. cutBack = .True.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator

6
src/mesh/FEM_utilities.f90
View File

@ -105,14 +105,14 @@ subroutine FEM_utilities_init
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
integer :: structOrder !< order of displacement shape functions integer :: structOrder !< order of displacement shape functions
character(len=pStringLen) :: & character(len=:), allocatable :: &
petsc_options petsc_options
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
structOrder = num_mesh%get_asInt('structOrder',defaultVal = 2) structOrder = num_mesh%get_asInt ('structOrder', defaultVal = 2)
petsc_options = num_mesh%get_asString('petsc_options', defaultVal='') petsc_options = num_mesh%get_asString('petsc_options', defaultVal='')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -134,7 +134,7 @@ subroutine FEM_utilities_init
&-mech_pc_type ml -mech_mg_levels_ksp_type chebyshev & &-mech_pc_type ml -mech_mg_levels_ksp_type chebyshev &
&-mech_mg_levels_pc_type sor -mech_pc_ml_nullspace user',ierr) &-mech_mg_levels_pc_type sor -mech_pc_ml_nullspace user',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr)

71
src/mesh/mesh_mech_FEM.f90
View File

@ -37,6 +37,18 @@ module mesh_mech_FEM
type(tSolutionParams) :: params type(tSolutionParams) :: params
type, private :: tNumerics
integer :: &
integrationOrder, & !< order of quadrature rule required
itmax
logical :: &
BBarStabilisation
real(pReal) :: &
eps_struct_atol, & !< absolute tolerance for mechanical equilibrium
eps_struct_rtol !< relative tolerance for mechanical equilibrium
end type tNumerics
type(tNumerics), private :: num
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
SNES :: mech_snes SNES :: mech_snes
@ -97,20 +109,21 @@ subroutine FEM_mech_init(fieldBC)
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
integer :: &
integrationOrder, & !< order of quadrature rule required
itmax
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'; flush(6)
!----------------------------------------------------------------------------- !-----------------------------------------------------------------------------
! read numerical parametes and do sanity checks ! read numerical parametes and do sanity checks
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) num%integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2)
num%itmax = num_mesh%get_asInt('itmax',defaultVal=250)
num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
num%eps_struct_atol = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal)
num%eps_struct_rtol = num_mesh%get_asFloat('eps_struct_rtol', defaultVal = 1.0e-4_pReal)
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
itmax = num_mesh%get_asInt('itmax',defaultVal=250) if (num%eps_struct_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_rtol')
if (itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%eps_struct_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_atol')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech mesh ! Setup FEM mech mesh
@ -119,9 +132,9 @@ subroutine FEM_mech_init(fieldBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech discretization ! Setup FEM mech discretization
qPoints = FEM_quadrature_points( dimPlex,integrationOrder)%p qPoints = FEM_quadrature_points( dimPlex,num%integrationOrder)%p
qWeights = FEM_quadrature_weights(dimPlex,integrationOrder)%p qWeights = FEM_quadrature_weights(dimPlex,num%integrationOrder)%p
nQuadrature = FEM_nQuadrature( dimPlex,integrationOrder) nQuadrature = FEM_nQuadrature( dimPlex,num%integrationOrder)
qPointsP => qPoints qPointsP => qPoints
qWeightsP => qWeights qWeightsP => qWeights
call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr) call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr)
@ -130,7 +143,7 @@ subroutine FEM_mech_init(fieldBC)
call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr) call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, & call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, &
integrationOrder,mechFE,ierr); CHKERRQ(ierr) num%integrationOrder,mechFE,ierr); CHKERRQ(ierr)
call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
nBasis = nBasis/nc nBasis = nBasis/nc
@ -221,7 +234,7 @@ subroutine FEM_mech_init(fieldBC)
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures
call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr) call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetTolerances(mech_snes,1.0,0.0,0.0,itmax,itmax,ierr) call SNESSetTolerances(mech_snes,1.0,0.0,0.0,num%itmax,num%itmax,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr)
@ -281,22 +294,10 @@ type(tSolutionState) function FEM_mech_solution( &
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
!--------------------------------------------------------------------------------------------------
integer :: &
itmax
class(tNode), pointer :: &
num_mesh
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
!--------------------------------------------------------------------------------------------------
! read numerical parameter and do sanity check
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
itmax = num_mesh%get_asInt('itmax',defaultVal=250)
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
!-------------------------------------------------------------------------------------------------
incInfo = incInfoIn incInfo = incInfoIn
FEM_mech_solution%converged =.false. FEM_mech_solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -311,7 +312,7 @@ type(tSolutionState) function FEM_mech_solution( &
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
FEM_mech_solution%converged = .false. FEM_mech_solution%converged = .false.
FEM_mech_solution%iterationsNeeded = itmax FEM_mech_solution%iterationsNeeded = num%itmax
else ! >= 1 proper convergence (or terminally ill) else ! >= 1 proper convergence (or terminally ill)
FEM_mech_solution%converged = .true. FEM_mech_solution%converged = .true.
call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr) call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr)
@ -349,12 +350,6 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
PetscInt :: bcSize PetscInt :: bcSize
IS :: bcPoints IS :: bcPoints
class(tNode), pointer :: &
num_mesh
logical :: BBarStabilisation
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2)) allocate(pcellJ(dimPlex**2))
@ -409,7 +404,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo enddo
if (BBarStabilisation) then if (num%BBarStabilisation) then
detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature))
do qPt = 1, nQuadrature do qPt = 1, nQuadrature
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = & materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = &
@ -498,12 +493,6 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
IS :: bcPoints IS :: bcPoints
class(tNode), pointer :: &
num_mesh
logical :: BBarStabilisation
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2)) allocate(pcellJ(dimPlex**2))
@ -560,7 +549,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
if (BBarStabilisation) then if (num%BBarStabilisation) then
F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex]) F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex])
FInv = math_inv33(F) FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
@ -577,7 +566,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
K_eA = K_eA + matmul(transpose(BMat),MatA) K_eA = K_eA + matmul(transpose(BMat),MatA)
endif endif
enddo enddo
if (BBarStabilisation) then if (num%BBarStabilisation) then
FInv = math_inv33(FAvg) FInv = math_inv33(FAvg)
K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + & K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + &
(matmul(matmul(transpose(BMatAvg), & (matmul(matmul(transpose(BMatAvg), &
@ -687,7 +676,7 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*1.0e-4_pReal,1.0e-10_pReal) divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*num%eps_struct_rtol,num%eps_struct_atol)
call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr) call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN