further cleaning of numerics.f90

This commit is contained in:
Sharan Roongta 2020-06-17 20:46:03 +02:00
parent d3f9e9f115
commit ac2539b305
9 changed files with 175 additions and 65 deletions

View File

@ -137,14 +137,23 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment
integer :: i, j, k, cell
integer :: i, j, k, cell, &
itmax !< maximum number of iterations
type(tSolutionState) :: solution
class(tNode), pointer :: &
num_generic
PetscInt :: devNull
PetscReal :: phi_min, phi_max, stagNorm, solnNorm
PetscErrorCode :: ierr
SNESConvergedReason :: reason
!-------------------------------------------------------------------
! reading numerical parameter and do sanity check
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmax = num_generic%get_asInt('itmax',defaultVal=250)
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
solution%converged =.false.
!--------------------------------------------------------------------------------------------------

View File

@ -410,23 +410,41 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
real(pReal) :: &
err_div, &
divTol, &
BCTol, &
eps_div_atol, &
eps_div_rtol, &
eps_stress_atol, &
eps_stress_rtol
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
num_grid, &
num_generic
!-----------------------------------------------------------------------------------
! 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)
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmin = num_generic%get_asInt('itmin',defaultVal=1)
itmax = num_generic%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
divTol = max(maxval(abs(P_av))*eps_div_rtol ,eps_div_atol)
BCTol = max(maxval(abs(P_av))*eps_stress_rtol,eps_stress_atol)
@ -474,7 +492,19 @@ subroutine formResidual(da_local,x_local, &
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal), dimension(3,3,3,3) :: devNull
integer :: &
itmin, &
itmax
class(tNode), pointer :: &
num_generic
!----------------------------------------------------------------------
! read numerical paramteters and do sanity checks
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmin = num_generic%get_asInt('itmin',defaultVal=1)
itmax = num_generic%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 SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr)

View File

@ -384,22 +384,39 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
real(pReal) :: &
divTol, &
BCTol, &
eps_div_atol, &
eps_div_rtol, &
eps_stress_atol, &
eps_stress_rtol
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
num_grid, &
num_generic
!-----------------------------------------------------------------------------------
! 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)
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmin = num_generic%get_asInt('itmin',defaultVal=1)
itmax = num_generic%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)
@ -445,6 +462,19 @@ subroutine formResidual(in, F, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
integer :: &
itmin, &
itmax
class(tNode), pointer :: &
num_generic
!----------------------------------------------------------------------
! read numerical paramteter and do sanity checks
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmin = num_generic%get_asInt('itmin',defaultVal=1)
itmax = num_generic%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 SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)

View File

@ -431,19 +431,25 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
real(pReal) :: &
curlTol, &
divTol, &
BCTol, &
eps_div_atol, &
eps_div_rtol, &
eps_curl_atol, &
eps_curl_rtol, &
eps_stress_atol, &
eps_stress_rtol
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
num_grid, &
num_generic
!-----------------------------------------------------------------------------------
! 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)
@ -452,6 +458,19 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
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)
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmin = num_generic%get_asInt('itmin',defaultVal=1)
itmax = num_generic%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)
@ -506,22 +525,30 @@ subroutine formResidual(in, FandF_tau, &
PetscObject :: dummy
PetscErrorCode :: ierr
class(tNode), pointer :: &
num_grid
num_grid, &
num_generic
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
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)
!-----------------------------------------------------------------------------------------------
! sanity checks
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmin = num_generic%get_asInt('itmin',defaultVal=1)
itmax = num_generic%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,&
XG_RANGE,YG_RANGE,ZG_RANGE)

View File

@ -11,9 +11,11 @@ module grid_thermal_spectral
use PETScsnes
use prec
use IO
use spectral_utilities
use discretization_grid
use thermal_conduction
use YAML_types
use numerics
use material
@ -132,14 +134,23 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment
integer :: i, j, k, cell
integer :: i, j, k, cell, &
itmax !< maximum number of iterations
type(tSolutionState) :: solution
class(tNode), pointer :: &
num_generic
PetscInt :: devNull
PetscReal :: T_min, T_max, stagNorm, solnNorm
PetscErrorCode :: ierr
SNESConvergedReason :: reason
!-------------------------------------------------------------------
! reading numerical parameter and do sanity check
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmax = num_generic%get_asInt('itmax',defaultVal=250)
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
solution%converged =.false.
!--------------------------------------------------------------------------------------------------

View File

@ -103,15 +103,15 @@ subroutine FEM_utilities_init
character(len=pStringLen) :: petsc_optionsOrder
class(tNode), pointer :: &
numerics_mesh
num_mesh
integer :: structOrder !< order of displacement shape functions
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
numerics_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
structOrder = numerics_mesh%get_asInt('structOrder',defaultVal = 2)
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
structOrder = num_mesh%get_asInt('structOrder',defaultVal = 2)
!--------------------------------------------------------------------------------------------------
! set debugging parameters

View File

@ -86,11 +86,11 @@ subroutine discretization_mesh_init(restart)
PetscErrorCode :: ierr
class(tNode), pointer :: &
numerics_mesh
num_mesh
integer :: integrationOrder !< order of quadrature rule required
numerics_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
integrationOrder = numerics_mesh%get_asInt('integrationorder',defaultVal = 2)
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2)
write(6,'(/,a)') ' <<<+- mesh init -+>>>'

View File

@ -96,13 +96,22 @@ subroutine FEM_mech_init(fieldBC)
PetscErrorCode :: ierr
class(tNode), pointer :: &
numerics_mesh
integer :: integrationOrder !< order of quadrature rule required
num_mesh, &
num_generic
integer :: &
integrationOrder, & !< order of quadrature rule required
itmax
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'; flush(6)
numerics_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
integrationOrder = numerics_mesh%get_asInt('integrationorder',defaultVal = 2)
!-----------------------------------------------------------------------------
! read numerical parametes and do sanity checks
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2)
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmax = num_generic%get_asInt('itmax',defaultVal=250)
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
!--------------------------------------------------------------------------------------------------
! Setup FEM mech mesh
@ -274,9 +283,21 @@ type(tSolutionState) function FEM_mech_solution( &
incInfoIn
!--------------------------------------------------------------------------------------------------
integer :: &
itmax
class(tNode), pointer :: &
num_generic
PetscErrorCode :: ierr
SNESConvergedReason :: reason
!--------------------------------------------------------------------------------------------------
! read numerical parameter and do sanity check
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmax = num_generic%get_asInt('itmax',defaultVal=250)
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
!-------------------------------------------------------------------------------------------------
incInfo = incInfoIn
FEM_mech_solution%converged =.false.
!--------------------------------------------------------------------------------------------------
@ -330,11 +351,11 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
IS :: bcPoints
class(tNode), pointer :: &
numerics_mesh
num_mesh
logical :: BBarStabilisation
numerics_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
BBarStabilisation = numerics_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2))
@ -479,11 +500,11 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
IS :: bcPoints
class(tNode), pointer :: &
numerics_mesh
num_mesh
logical :: BBarStabilisation
numerics_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
BBarStabilisation = numerics_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
num_mesh => numerics_root%get('mesh',defaultVal=emptyDict)
BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2))

View File

@ -25,25 +25,19 @@ module numerics
worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only)
integer(4), protected, public :: &
DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
!--------------------------------------------------------------------------------------------------
! field parameters:
integer, protected, public :: &
itmax = 250, & !< maximum number of iterations
itmin = 1 !< minimum number of iterations
!--------------------------------------------------------------------------------------------------
! spectral parameters:
#ifdef Grid
character(len=pStringLen), protected, public :: &
petsc_options = ''
petsc_options
#endif
!--------------------------------------------------------------------------------------------------
! Mesh parameters:
#ifdef Mesh
character(len=pStringLen), protected, public :: &
petsc_options = ''
petsc_options
#endif
public :: numerics_init
@ -58,7 +52,7 @@ contains
subroutine numerics_init
!$ integer :: gotDAMASK_NUM_THREADS = 1
integer :: i, ierr
integer :: ierr
character(len=:), allocatable :: &
numerics_input, &
numerics_inFlow
@ -95,8 +89,6 @@ subroutine numerics_init
!--------------------------------------------------------------------------------------------------
! grid parameters
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
itmax = num_grid%get_asInt('itmax',defaultVal=250)
itmin = num_grid%get_asInt('itmin',defaultVal=1)
#ifdef PETSC
petsc_options = num_grid%get_asString('petsc_options',defaultVal = '')
#endif
@ -109,21 +101,11 @@ subroutine numerics_init
! openMP parameter
!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
!--------------------------------------------------------------------------------------------------
! field parameters
write(6,'(a24,1x,i8)') ' itmax: ',itmax
write(6,'(a24,1x,i8)') ' itmin: ',itmin
!--------------------------------------------------------------------------------------------------
#ifdef PETSC
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options)
#endif
!--------------------------------------------------------------------------------------------------
! sanity checks
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
end subroutine numerics_init
end module numerics