diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 9113fa12f..2454abaab 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -83,7 +83,7 @@ module crystallite iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp nState, & !< state loop limit nStress !< stress loop limit - character(len=pStringLen) :: & + character(len=:), allocatable :: & integrator !< integrator scheme real(pReal) :: & 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') - select case(trim(num%integrator)) + select case(num%integrator) case('FPI') integrateState => integrateStateFPI case('Euler') diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index bdd76b423..775bc913d 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -30,6 +30,23 @@ module grid_mech_FEM ! derived types 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 DM, private :: mech_grid @@ -97,8 +114,7 @@ subroutine grid_mech_FEM_init PetscErrorCode :: ierr integer(HID_T) :: fileHandle, groupHandle character(len=pStringLen) :: & - fileName, & - petsc_options + fileName class(tNode), pointer :: & num_grid 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) !------------------------------------------------------------------------------------------------- -! read numerical parameter +! read numerical parameter and do sanity checks 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 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) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -419,49 +449,23 @@ 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, & !< 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 + BCTol -!----------------------------------------------------------------------------------- -! 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 - divTol = max(maxval(abs(P_av))*eps_div_rtol ,eps_div_atol) - BCTol = max(maxval(abs(P_av))*eps_stress_rtol,eps_stress_atol) + 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 >= itmin .and. & + if ((totalIter >= num%itmin .and. & all([ err_div/divTol, & err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 - elseif (totalIter >= itmax) then + elseif (totalIter >= num%itmax) then reason = -1 else reason = 0 @@ -499,19 +503,6 @@ 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_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 SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) @@ -522,7 +513,7 @@ subroutine formResidual(da_local,x_local, & ! begin of new iteration newIteration: if (totalIter <= PETScIter) then 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) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 61ac18779..b1c4d085a 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -32,9 +32,19 @@ module grid_mech_spectral_basic type, private :: tNumerics 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 - type(tNumerics) :: num ! numerics parameters. Better name? + type(tNumerics), private :: num ! numerics parameters. Better name? !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -98,8 +108,7 @@ subroutine grid_mech_spectral_basic_init integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: & - fileName, & - petsc_options + fileName 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' !------------------------------------------------------------------------------------------------- -! read numerical parameters +! read numerical parameters and do sanity checks 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 call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -389,47 +413,19 @@ 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, & !< 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 + BCTol -!----------------------------------------------------------------------------------- -! 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) + 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) - 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') -!------------------------------------------------------------------------------------ - - 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. & + if ((totalIter >= num%itmin .and. & all([ err_div/divTol, & err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 - elseif (totalIter >= itmax) then + elseif (totalIter >= num%itmax) then reason = -1 else reason = 0 @@ -466,19 +462,6 @@ subroutine formResidual(in, F, & nfuncs PetscObject :: dummy 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 SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) @@ -488,7 +471,7 @@ subroutine formResidual(in, F, & ! begin of new iteration newIteration: if (totalIter <= PETScIter) then 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) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index a32216615..87cb69579 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -33,9 +33,24 @@ module grid_mech_spectral_polarisation type, private :: tNumerics 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 - type(tNumerics) :: num ! numerics parameters. Better name? + type(tNumerics), private :: num ! numerics parameters. Better name? !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -107,8 +122,7 @@ subroutine grid_mech_spectral_polarisation_init integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: & - fileName, & - petsc_options + fileName write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(6) @@ -118,14 +132,36 @@ subroutine grid_mech_spectral_polarisation_init !------------------------------------------------------------------------------------------------- ! read numerical parameters 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 call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -436,56 +472,22 @@ 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, & !< 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 + BCTol -!----------------------------------------------------------------------------------- -! 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_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) + curlTol = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol ,num%eps_curl_atol) + 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) - 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_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. & + if ((totalIter >= num%itmin .and. & all([ err_div /divTol, & err_curl/curlTol, & err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 - elseif (totalIter >= itmax) then + elseif (totalIter >= num%itmax) then reason = -1 else reason = 0 @@ -527,27 +529,9 @@ subroutine formResidual(in, FandF_tau, & nfuncs PetscObject :: dummy PetscErrorCode :: ierr - class(tNode), pointer :: & - num_grid - 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 + integer :: & + i, j, k, e -!-------------------------------------------------------------------------------------------------- -! 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,& @@ -570,7 +554,7 @@ subroutine formResidual(in, FandF_tau, & ! begin of new iteration newIteration: if (totalIter <= PETScIter) then 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) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' 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 do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) tensorField_real(1:3,1:3,i,j,k) = & - 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%polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + 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)) enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space 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 !-------------------------------------------------------------------------------------------------- ! 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 call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & - F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) + F - residual_F_tau/num%polarBeta,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index cc9d3dd65..08e5d633c 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -121,10 +121,10 @@ module spectral_utilities character(len=:), allocatable :: & spectral_derivative, & !< approximation used for derivatives in Fourier space FFTW_plan_mode, & !< FFTW plan mode, see www.fftw.org - PETSc_options + petsc_options end type tNumerics - type(tNumerics) :: num ! numerics parameters. Better name? + type(tNumerics), private :: num ! numerics parameters. Better name? enum, bind(c); enumerator :: & DERIVATIVE_CONTINUOUS_ID, & @@ -189,8 +189,6 @@ subroutine spectral_utilities_init scalarSize = 1_C_INTPTR_T, & vecSize = 3_C_INTPTR_T, & tensorSize = 9_C_INTPTR_T - character(len=pStringLen) :: & - petsc_options class(tNode), pointer :: & num_grid @@ -220,13 +218,13 @@ subroutine spectral_utilities_init ' add more using the PETSc_Options keyword in numerics.config '; flush(6) 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) CHKERRQ(ierr) if(debugPETSc) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) grid1Red = grid(1)/2 + 1 diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index 8945c2c18..a4a187c38 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -34,8 +34,8 @@ module discretization_marc mesh_unitlength !< physical length of one unit in mesh integer, dimension(:), allocatable, public :: & - mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID - mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID + mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID + mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID public :: & discretization_marc_init diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 72fa53566..185f4c375 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -49,7 +49,6 @@ program DAMASK_mesh i, & errorID, & 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 currentLoadcase = 0, & !< current load case currentFace = 0, & @@ -57,7 +56,6 @@ program DAMASK_mesh totalIncsCounter = 0, & !< total # of increments statUnit = 0, & !< file unit for statistics output stagIter, & - stagItMax, & !< max number of field level staggered iterations component class(tNode), pointer :: & num_mesh @@ -65,6 +63,14 @@ program DAMASK_mesh character(len=pStringLen) :: & incInfo, & 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(tSolutionState), allocatable, dimension(:) :: solres PetscInt :: faceSet, currentFaceSet, field, dimPlex @@ -81,11 +87,11 @@ program DAMASK_mesh !--------------------------------------------------------------------- ! reading field information from numerics file and do sanity checks num_mesh => numerics_root%get('mesh', defaultVal=emptyDict) - stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10) - maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3) + num%stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10) + num%maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3) - if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') - if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') + if (num%stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') + 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 call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRA(ierr) !< dimension of mesh (2D or 3D) @@ -327,7 +333,7 @@ program DAMASK_mesh enddo stagIter = stagIter + 1 - stagIterate = stagIter < stagItMax & + stagIterate = stagIter < num%stagItMax & .and. all(solres(:)%converged) & .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration enddo @@ -335,7 +341,7 @@ program DAMASK_mesh ! check solution cutBack = .False. 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' cutBack = .True. stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index a8479913e..a6236175c 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -105,14 +105,14 @@ subroutine FEM_utilities_init class(tNode), pointer :: & num_mesh integer :: structOrder !< order of displacement shape functions - character(len=pStringLen) :: & + character(len=:), allocatable :: & petsc_options PetscErrorCode :: ierr write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>' 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='') !-------------------------------------------------------------------------------------------------- @@ -134,7 +134,7 @@ subroutine FEM_utilities_init &-mech_pc_type ml -mech_mg_levels_ksp_type chebyshev & &-mech_mg_levels_pc_type sor -mech_pc_ml_nullspace user',ierr) CHKERRQ(ierr) - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,ierr) CHKERRQ(ierr) write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 715b02c0d..918914432 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -37,6 +37,18 @@ module mesh_mech_FEM 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 SNES :: mech_snes @@ -97,20 +109,21 @@ subroutine FEM_mech_init(fieldBC) class(tNode), pointer :: & num_mesh - integer :: & - integrationOrder, & !< order of quadrature rule required - itmax write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'; flush(6) !----------------------------------------------------------------------------- ! read numerical parametes and do sanity checks 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) - itmax = num_mesh%get_asInt('itmax',defaultVal=250) - if (itmax <= 1) call IO_error(301,ext_msg='itmax') + if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') + if (num%eps_struct_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_rtol') + if (num%eps_struct_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_atol') !-------------------------------------------------------------------------------------------------- ! Setup FEM mech mesh @@ -119,9 +132,9 @@ subroutine FEM_mech_init(fieldBC) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech discretization - qPoints = FEM_quadrature_points( dimPlex,integrationOrder)%p - qWeights = FEM_quadrature_weights(dimPlex,integrationOrder)%p - nQuadrature = FEM_nQuadrature( dimPlex,integrationOrder) + qPoints = FEM_quadrature_points( dimPlex,num%integrationOrder)%p + qWeights = FEM_quadrature_weights(dimPlex,num%integrationOrder)%p + nQuadrature = FEM_nQuadrature( dimPlex,num%integrationOrder) qPointsP => qPoints qWeightsP => qWeights 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) CHKERRQ(ierr) 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 PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) 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 SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,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) call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) @@ -281,22 +294,10 @@ type(tSolutionState) function FEM_mech_solution( & character(len=*), intent(in) :: & incInfoIn -!-------------------------------------------------------------------------------------------------- - integer :: & - itmax - class(tNode), pointer :: & - num_mesh PetscErrorCode :: ierr 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 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 FEM_mech_solution%converged = .false. - FEM_mech_solution%iterationsNeeded = itmax + FEM_mech_solution%iterationsNeeded = num%itmax else ! >= 1 proper convergence (or terminally ill) FEM_mech_solution%converged = .true. 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 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(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) = & reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) 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)) do qPt = 1, nQuadrature 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 - 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(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), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & 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]) FInv = math_inv33(F) 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) endif enddo - if (BBarStabilisation) then + if (num%BBarStabilisation) then FInv = math_inv33(FAvg) K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + & (matmul(matmul(transpose(BMatAvg), & @@ -687,7 +676,7 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm !-------------------------------------------------------------------------------------------------- ! 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) CHKERRQ(ierr) if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN