2011-04-07 12:50:28 +05:30
! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
2011-04-04 19:39:54 +05:30
!
! This file is part of DAMASK,
2011-04-07 12:50:28 +05:30
! the Düsseldorf Advanced MAterial Simulation Kit.
2011-04-04 19:39:54 +05:30
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Managing of parameters related to numerics
!--------------------------------------------------------------------------------------------------
2012-03-07 15:37:29 +05:30
module numerics
2013-01-10 03:49:32 +05:30
use prec , only : &
pInt , &
pReal
2009-06-15 18:41:21 +05:30
2013-01-10 03:49:32 +05:30
implicit none
private
character ( len = 64 ) , parameter , private :: &
numerics_configFile = 'numerics.config' !< name of configuration file
2009-06-15 18:41:21 +05:30
2013-01-10 03:49:32 +05:30
integer ( pInt ) , protected , public :: &
2012-10-02 20:56:56 +05:30
iJacoStiffness = 1_pInt , & !< frequency of stiffness update
2012-09-13 15:18:38 +05:30
iJacoLpresiduum = 1_pInt , & !< frequency of Jacobian update of residuum in Lp
nHomog = 20_pInt , & !< homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog")
nMPstate = 10_pInt , & !< materialpoint state loop limit
nCryst = 20_pInt , & !< crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")
nState = 10_pInt , & !< state loop limit
nStress = 40_pInt , & !< stress loop limit
2012-10-02 20:56:56 +05:30
pert_method = 1_pInt !< method used in perturbation technique for tangent
2013-01-10 03:49:32 +05:30
integer ( pInt ) , public :: numerics_integrationMode = 0_pInt !< integrationMode 1 = central solution ; integrationMode 2 = perturbation, Default 0: undefined, is not read from file
integer ( pInt ) , dimension ( 2 ) , protected , public :: &
2012-10-02 20:56:56 +05:30
numerics_integrator = 1_pInt !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states
2013-01-10 03:49:32 +05:30
real ( pReal ) , protected , public :: &
2012-10-02 20:56:56 +05:30
relevantStrain = 1.0e-7_pReal , & !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
2012-09-13 15:18:38 +05:30
defgradTolerance = 1.0e-7_pReal , & !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
pert_Fg = 1.0e-7_pReal , & !< strain perturbation for FEM Jacobi
subStepMinCryst = 1.0e-3_pReal , & !< minimum (relative) size of sub-step allowed during cutback in crystallite
subStepMinHomog = 1.0e-3_pReal , & !< minimum (relative) size of sub-step allowed during cutback in homogenization
subStepSizeCryst = 0.25_pReal , & !< size of first substep when cutback in crystallite
subStepSizeHomog = 0.25_pReal , & !< size of first substep when cutback in homogenization
stepIncreaseCryst = 1.5_pReal , & !< increase of next substep size when previous substep converged in crystallite
stepIncreaseHomog = 1.5_pReal , & !< increase of next substep size when previous substep converged in homogenization
rTol_crystalliteState = 1.0e-6_pReal , & !< relative tolerance in crystallite state loop
rTol_crystalliteTemperature = 1.0e-6_pReal , & !< relative tolerance in crystallite temperature loop
rTol_crystalliteStress = 1.0e-6_pReal , & !< relative tolerance in crystallite stress loop
aTol_crystalliteStress = 1.0e-8_pReal , & !< absolute tolerance in crystallite stress loop, Default 1.0e-8: residuum is in Lp and hence strain is on this order
2012-09-24 21:52:25 +05:30
numerics_unitlength = 1.0_pReal , & !< determines the physical length of one computational length unit
2012-09-13 15:18:38 +05:30
absTol_RGC = 1.0e+4_pReal , & !< absolute tolerance of RGC residuum
relTol_RGC = 1.0e-3_pReal , & !< relative tolerance of RGC residuum
absMax_RGC = 1.0e+10_pReal , & !< absolute maximum of RGC residuum
relMax_RGC = 1.0e+2_pReal , & !< relative maximum of RGC residuum
pPert_RGC = 1.0e-7_pReal , & !< perturbation for computing RGC penalty tangent
xSmoo_RGC = 1.0e-5_pReal , & !< RGC penalty smoothing parameter (hyperbolic tangent)
viscPower_RGC = 1.0e+0_pReal , & !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model)
viscModus_RGC = 0.0e+0_pReal , & !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied
refRelaxRate_RGC = 1.0e-3_pReal , & !< reference relaxation rate in RGC viscosity
maxdRelax_RGC = 1.0e+0_pReal , & !< threshold of maximum relaxation vector increment (if exceed this then cutback)
maxVolDiscr_RGC = 1.0e-5_pReal , & !< threshold of maximum volume discrepancy allowed
volDiscrMod_RGC = 1.0e+12_pReal , & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
volDiscrPow_RGC = 5.0_pReal !< powerlaw penalty for volume discrepancy
2013-01-10 03:49:32 +05:30
logical , protected , public :: &
2012-11-28 00:06:55 +05:30
analyticJaco = . false . , & !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations
numerics_timeSyncing = . false . !< flag indicating if time synchronization in crystallite is used for nonlocal plasticity
2012-06-15 21:40:21 +05:30
!* Random seeding parameters
2013-01-10 03:49:32 +05:30
integer ( pInt ) , protected , public :: &
2012-10-02 20:56:56 +05:30
fixedSeed = 0_pInt !< fixed seeding for pseudo-random number generator, Default 0: use random seed
2012-06-15 21:40:21 +05:30
!* OpenMP variable
2013-01-10 03:49:32 +05:30
integer ( pInt ) , protected , public :: &
2012-10-02 20:56:56 +05:30
DAMASK_NumThreadsInt = 0_pInt !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
2012-06-15 21:40:21 +05:30
2010-10-13 21:34:44 +05:30
!* spectral parameters:
2012-06-15 21:40:21 +05:30
#ifdef Spectral
2013-01-10 03:49:32 +05:30
real ( pReal ) , protected , public :: &
2013-03-13 00:18:28 +05:30
err_div_tol = 1.0e-5_pReal , & !< Div(P)/avg(P)*meter
2012-09-13 15:18:38 +05:30
err_stress_tolrel = 0.01_pReal , & !< relative tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress
err_stress_tolabs = huge ( 1.0_pReal ) , & !< absolute tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress
2013-03-13 00:18:28 +05:30
err_f_tol = 1.0e-6_pReal , &
err_p_tol = 1.0e-5_pReal , &
2012-09-13 15:18:38 +05:30
fftw_timelimit = - 1.0_pReal , & !< sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit
2013-02-28 23:05:02 +05:30
rotation_tol = 1.0e-12_pReal , & !< tolerance of rotation specified in loadcase, Default 1.0e-12: first guess
polarAlpha = 1.0_pReal , & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
polarBeta = 1.0_pReal !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
2013-01-10 03:49:32 +05:30
character ( len = 64 ) , private :: &
2012-10-02 20:56:56 +05:30
fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag
2013-01-10 03:49:32 +05:30
character ( len = 64 ) , protected , public :: &
2012-09-13 15:18:38 +05:30
myspectralsolver = 'basic' , & !< spectral solution method
myfilter = 'none' !< spectral filtering method
2013-01-10 03:49:32 +05:30
character ( len = 1024 ) , protected , public :: &
2012-12-15 23:37:49 +05:30
petsc_options = ' - snes_type ngmres &
2013-01-03 21:47:23 +05:30
& - snes_ngmres_anderson '
2013-01-10 03:49:32 +05:30
integer ( pInt ) , protected , public :: &
2012-10-02 20:56:56 +05:30
fftw_planner_flag = 32_pInt , & !< conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw
2013-03-15 04:40:02 +05:30
itmax = 250_pInt , & !< maximum number of iterations
2012-09-13 15:18:38 +05:30
itmin = 2_pInt , & !< minimum number of iterations
2012-12-14 20:48:04 +05:30
maxCutBack = 3_pInt , & !< max number of cut backs
2013-01-16 16:10:53 +05:30
regridMode = 0_pInt , & !< 0: no regrid; 1: regrid if DAMASK doesn't converge; 2: regrid if DAMASK or BVP Solver doesn't converge
2013-03-13 00:18:28 +05:30
divergence_correction = 2_pInt !< correct divergence calculation in fourier space 0: no correction, 1: dimension scaled to 1, 2: dimension scaled to Npoints, 3: dimension scaled to sqrt(Npoints)
2013-01-10 03:49:32 +05:30
logical , protected , public :: &
2012-10-02 20:56:56 +05:30
memory_efficient = . true . , & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate
2012-09-13 15:18:38 +05:30
update_gamma = . false . !< update gamma operator with current stiffness, Default .false.: use initial stiffness
2009-06-15 18:41:21 +05:30
2012-10-02 20:56:56 +05:30
#endif
2011-02-07 20:05:42 +05:30
2013-01-10 03:49:32 +05:30
public :: numerics_init
2012-10-02 20:56:56 +05:30
contains
2009-11-10 19:06:27 +05:30
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief reads in parameters from numerics.config and sets openMP related parameters. Also does
! a sanity check
!--------------------------------------------------------------------------------------------------
2012-03-07 15:37:29 +05:30
subroutine numerics_init
2012-06-12 15:14:05 +05:30
use , intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
2012-12-15 23:37:49 +05:30
use IO , only : &
IO_error , &
IO_open_file_stat , &
IO_isBlank , &
IO_stringPos , &
IO_stringValue , &
IO_lc , &
IO_floatValue , &
IO_intValue , &
2013-02-25 22:04:59 +05:30
IO_warning , &
IO_timeStamp
2012-12-15 23:37:49 +05:30
2012-12-15 21:51:10 +05:30
#ifndef Marc
!$ use OMP_LIB, only: omp_set_num_threads ! Use the standard conforming module file for omp if not using MSC.Marc
#endif
2012-06-15 21:40:21 +05:30
implicit none
2012-12-15 21:51:10 +05:30
#ifdef Marc
2013-02-11 15:14:17 +05:30
!$ include "omp_lib.h" ! use the not F90 standard conforming include file to prevent crashes with some versions of MSC.Marc
2012-12-15 21:51:10 +05:30
#endif
2012-03-07 15:37:29 +05:30
integer ( pInt ) , parameter :: fileunit = 300_pInt , &
maxNchunks = 2_pInt
!$ integer :: gotDAMASK_NUM_THREADS = 1
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: positions
character ( len = 64 ) :: tag
character ( len = 1024 ) :: line
2012-06-12 15:14:05 +05:30
!$ character(len=6) DAMASK_NumThreadsString ! environment variable DAMASK_NUM_THREADS
2009-08-27 21:00:40 +05:30
2013-02-11 15:14:17 +05:30
write ( 6 , '(/,a)' ) ' <<<+- numerics init -+>>>'
write ( 6 , '(a)' ) ' $Id$'
2013-02-25 22:04:59 +05:30
write ( 6 , '(a16,a)' ) ' Current time : ' , IO_timeStamp ( )
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
2012-06-15 21:40:21 +05:30
2012-01-31 01:46:19 +05:30
!$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS...
2013-02-21 03:26:59 +05:30
!$ if(gotDAMASK_NUM_THREADS /= 0) &
!$ call IO_warning(35_pInt,ext_msg='BEGIN:'//DAMASK_NumThreadsString//':END')
2012-01-31 01:46:19 +05:30
!$ read(DAMASK_NumThreadsString,'(i6)') DAMASK_NumThreadsInt ! ...convert it to integer...
2012-02-16 00:28:38 +05:30
!$ if (DAMASK_NumThreadsInt < 1_pInt) DAMASK_NumThreadsInt = 1_pInt ! ...ensure that its at least one...
2012-01-30 19:22:41 +05:30
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! ...and use it as number of threads for parallel execution
2010-12-02 16:34:29 +05:30
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! try to open the config file
fileExists : if ( IO_open_file_stat ( fileunit , numerics_configFile ) ) then
write ( 6 , '(a,/)' ) ' using values from config file'
2012-02-13 23:11:27 +05:30
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! read variables from config file and overwrite default parameters if keyword is present
2012-03-07 15:37:29 +05:30
line = ''
do
read ( fileunit , '(a1024)' , END = 100 ) line
2013-02-11 15:14:17 +05:30
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
2012-03-07 15:37:29 +05:30
positions = IO_stringPos ( line , maxNchunks )
2013-02-11 15:14:17 +05:30
tag = IO_lc ( IO_stringValue ( line , positions , 1_pInt ) ) ! extract key
2012-03-07 15:37:29 +05:30
select case ( tag )
case ( 'relevantstrain' )
relevantStrain = IO_floatValue ( line , positions , 2_pInt )
case ( 'defgradtolerance' )
defgradTolerance = IO_floatValue ( line , positions , 2_pInt )
case ( 'ijacostiffness' )
iJacoStiffness = IO_intValue ( line , positions , 2_pInt )
case ( 'ijacolpresiduum' )
iJacoLpresiduum = IO_intValue ( line , positions , 2_pInt )
case ( 'pert_fg' )
pert_Fg = IO_floatValue ( line , positions , 2_pInt )
case ( 'pert_method' )
pert_method = IO_intValue ( line , positions , 2_pInt )
case ( 'nhomog' )
nHomog = IO_intValue ( line , positions , 2_pInt )
case ( 'nmpstate' )
nMPstate = IO_intValue ( line , positions , 2_pInt )
case ( 'ncryst' )
nCryst = IO_intValue ( line , positions , 2_pInt )
case ( 'nstate' )
nState = IO_intValue ( line , positions , 2_pInt )
case ( 'nstress' )
nStress = IO_intValue ( line , positions , 2_pInt )
case ( 'substepmincryst' )
subStepMinCryst = IO_floatValue ( line , positions , 2_pInt )
case ( 'substepsizecryst' )
subStepSizeCryst = IO_floatValue ( line , positions , 2_pInt )
case ( 'stepincreasecryst' )
stepIncreaseCryst = IO_floatValue ( line , positions , 2_pInt )
case ( 'substepminhomog' )
subStepMinHomog = IO_floatValue ( line , positions , 2_pInt )
case ( 'substepsizehomog' )
subStepSizeHomog = IO_floatValue ( line , positions , 2_pInt )
case ( 'stepincreasehomog' )
stepIncreaseHomog = IO_floatValue ( line , positions , 2_pInt )
case ( 'rtol_crystallitestate' )
rTol_crystalliteState = IO_floatValue ( line , positions , 2_pInt )
case ( 'rtol_crystallitetemperature' )
rTol_crystalliteTemperature = IO_floatValue ( line , positions , 2_pInt )
case ( 'rtol_crystallitestress' )
rTol_crystalliteStress = IO_floatValue ( line , positions , 2_pInt )
case ( 'atol_crystallitestress' )
aTol_crystalliteStress = IO_floatValue ( line , positions , 2_pInt )
case ( 'integrator' )
numerics_integrator ( 1 ) = IO_intValue ( line , positions , 2_pInt )
case ( 'integratorstiffness' )
numerics_integrator ( 2 ) = IO_intValue ( line , positions , 2_pInt )
case ( 'analyticjaco' )
analyticJaco = IO_intValue ( line , positions , 2_pInt ) > 0_pInt
2012-11-28 00:06:55 +05:30
case ( 'timesyncing' )
numerics_timeSyncing = IO_intValue ( line , positions , 2_pInt ) > 0_pInt
2012-09-24 21:52:25 +05:30
case ( 'unitlength' )
numerics_unitlength = IO_floatValue ( line , positions , 2_pInt )
2009-07-31 17:32:20 +05:30
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! RGC parameters
2012-03-07 15:37:29 +05:30
case ( 'atol_rgc' )
absTol_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'rtol_rgc' )
relTol_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'amax_rgc' )
absMax_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'rmax_rgc' )
relMax_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'perturbpenalty_rgc' )
pPert_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'relevantmismatch_rgc' )
xSmoo_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'viscositypower_rgc' )
viscPower_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'viscositymodulus_rgc' )
viscModus_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'refrelaxationrate_rgc' )
refRelaxRate_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'maxrelaxation_rgc' )
maxdRelax_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'maxvoldiscrepancy_rgc' )
maxVolDiscr_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'voldiscrepancymod_rgc' )
volDiscrMod_RGC = IO_floatValue ( line , positions , 2_pInt )
case ( 'discrepancypower_rgc' )
volDiscrPow_RGC = IO_floatValue ( line , positions , 2_pInt )
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! random seeding parameters
2012-06-15 21:40:21 +05:30
case ( 'fixed_seed' )
fixedSeed = IO_intValue ( line , positions , 2_pInt )
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! spectral parameters
2012-06-15 21:40:21 +05:30
#ifdef Spectral
2012-03-07 15:37:29 +05:30
case ( 'err_div_tol' )
err_div_tol = IO_floatValue ( line , positions , 2_pInt )
case ( 'err_stress_tolrel' )
err_stress_tolrel = IO_floatValue ( line , positions , 2_pInt )
2012-04-11 18:27:25 +05:30
case ( 'err_stress_tolabs' )
err_stress_tolabs = IO_floatValue ( line , positions , 2_pInt )
2012-03-07 15:37:29 +05:30
case ( 'itmax' )
itmax = IO_intValue ( line , positions , 2_pInt )
case ( 'itmin' )
itmin = IO_intValue ( line , positions , 2_pInt )
2012-09-13 15:18:38 +05:30
case ( 'maxcutback' )
maxCutBack = IO_intValue ( line , positions , 2_pInt )
2012-12-14 20:48:04 +05:30
case ( 'regridmode' )
regridMode = IO_intValue ( line , positions , 2_pInt )
2012-03-07 15:37:29 +05:30
case ( 'memory_efficient' )
memory_efficient = IO_intValue ( line , positions , 2_pInt ) > 0_pInt
case ( 'fftw_timelimit' )
fftw_timelimit = IO_floatValue ( line , positions , 2_pInt )
case ( 'fftw_plan_mode' )
2012-10-02 20:56:56 +05:30
fftw_plan_mode = IO_lc ( IO_stringValue ( line , positions , 2_pInt ) )
2012-08-07 22:53:13 +05:30
case ( 'myfilter' )
2012-10-02 20:56:56 +05:30
myfilter = IO_lc ( IO_stringValue ( line , positions , 2_pInt ) )
2012-03-07 15:37:29 +05:30
case ( 'rotation_tol' )
rotation_tol = IO_floatValue ( line , positions , 2_pInt )
case ( 'divergence_correction' )
2013-01-16 16:10:53 +05:30
divergence_correction = IO_intValue ( line , positions , 2_pInt )
2012-03-07 15:37:29 +05:30
case ( 'update_gamma' )
update_gamma = IO_intValue ( line , positions , 2_pInt ) > 0_pInt
2012-08-29 00:49:47 +05:30
#ifdef PETSc
case ( 'petsc_options' )
petsc_options = trim ( line ( positions ( 4 ) : ) )
case ( 'myspectralsolver' )
2012-10-02 20:56:56 +05:30
myspectralsolver = IO_lc ( IO_stringValue ( line , positions , 2_pInt ) )
2012-08-29 00:49:47 +05:30
case ( 'err_f_tol' )
err_f_tol = IO_floatValue ( line , positions , 2_pInt )
case ( 'err_p_tol' )
err_p_tol = IO_floatValue ( line , positions , 2_pInt )
2013-02-28 23:05:02 +05:30
case ( 'polaralpha' )
polarAlpha = IO_floatValue ( line , positions , 2_pInt )
case ( 'polarbeta' )
polarBeta = IO_floatValue ( line , positions , 2_pInt )
2012-08-29 00:49:47 +05:30
#endif
#ifndef PETSc
2013-02-28 23:05:02 +05:30
case ( 'myspectralsolver' , 'petsc_options' , 'err_f_tol' , 'err_p_tol' , &
'polaralpha' , 'polarBeta' ) ! found PETSc parameter, but compiled without PETSc
2012-08-29 00:49:47 +05:30
call IO_warning ( 41_pInt , ext_msg = tag )
#endif
2012-06-15 21:40:21 +05:30
#endif
#ifndef Spectral
2013-02-11 15:14:17 +05:30
case ( 'err_div_tol' , 'err_stress_tolrel' , 'err_stress_tolabs' , & ! found spectral parameter for FEM build
2012-07-25 19:31:39 +05:30
'itmax' , 'itmin' , 'memory_efficient' , 'fftw_timelimit' , 'fftw_plan_mode' , 'myspectralsolver' , &
2012-08-29 00:49:47 +05:30
'rotation_tol' , 'divergence_correction' , 'update_gamma' , 'petsc_options' , 'myfilter' , &
2013-02-28 23:05:02 +05:30
'err_f_tol' , 'err_p_tol' , 'maxcutback' , 'polaralpha' , 'polarbeta' )
2012-06-15 21:40:21 +05:30
call IO_warning ( 40_pInt , ext_msg = tag )
#endif
2013-02-11 15:14:17 +05:30
case default ! found unknown keyword
2012-03-07 15:37:29 +05:30
call IO_error ( 300_pInt , ext_msg = tag )
endselect
enddo
100 close ( fileunit )
2013-02-11 15:14:17 +05:30
else fileExists
write ( 6 , '(a,/)' ) ' using standard values'
endif fileExists
2012-06-15 21:40:21 +05:30
#ifdef Spectral
2013-02-11 15:14:17 +05:30
select case ( IO_lc ( fftw_plan_mode ) ) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
case ( 'estimate' , 'fftw_estimate' ) ! ordered from slow execution (but fast plan creation) to fast execution
2012-03-07 15:37:29 +05:30
fftw_planner_flag = 64_pInt
case ( 'measure' , 'fftw_measure' )
fftw_planner_flag = 0_pInt
case ( 'patient' , 'fftw_patient' )
fftw_planner_flag = 32_pInt
case ( 'exhaustive' , 'fftw_exhaustive' )
fftw_planner_flag = 8_pInt
case default
call IO_warning ( warning_ID = 47_pInt , ext_msg = trim ( IO_lc ( fftw_plan_mode ) ) )
fftw_planner_flag = 32_pInt
end select
2012-06-15 21:40:21 +05:30
#endif
2012-11-28 00:06:55 +05:30
2013-02-11 15:14:17 +05:30
numerics_timeSyncing = numerics_timeSyncing . and . all ( numerics_integrator == 2_pInt ) ! timeSyncing only allowed for explicit Euler integrator
2012-11-28 00:06:55 +05:30
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! writing parameters to output file
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,es8.1)' ) ' relevantStrain: ' , relevantStrain
write ( 6 , '(a24,1x,es8.1)' ) ' defgradTolerance: ' , defgradTolerance
write ( 6 , '(a24,1x,i8)' ) ' iJacoStiffness: ' , iJacoStiffness
write ( 6 , '(a24,1x,i8)' ) ' iJacoLpresiduum: ' , iJacoLpresiduum
write ( 6 , '(a24,1x,es8.1)' ) ' pert_Fg: ' , pert_Fg
write ( 6 , '(a24,1x,i8)' ) ' pert_method: ' , pert_method
write ( 6 , '(a24,1x,i8)' ) ' nCryst: ' , nCryst
write ( 6 , '(a24,1x,es8.1)' ) ' subStepMinCryst: ' , subStepMinCryst
write ( 6 , '(a24,1x,es8.1)' ) ' subStepSizeCryst: ' , subStepSizeCryst
write ( 6 , '(a24,1x,es8.1)' ) ' stepIncreaseCryst: ' , stepIncreaseCryst
write ( 6 , '(a24,1x,i8)' ) ' nState: ' , nState
write ( 6 , '(a24,1x,i8)' ) ' nStress: ' , nStress
write ( 6 , '(a24,1x,es8.1)' ) ' rTol_crystalliteState: ' , rTol_crystalliteState
write ( 6 , '(a24,1x,es8.1)' ) ' rTol_crystalliteTemp: ' , rTol_crystalliteTemperature
write ( 6 , '(a24,1x,es8.1)' ) ' rTol_crystalliteStress: ' , rTol_crystalliteStress
write ( 6 , '(a24,1x,es8.1)' ) ' aTol_crystalliteStress: ' , aTol_crystalliteStress
write ( 6 , '(a24,2(1x,i8))' ) ' integrator: ' , numerics_integrator
2012-11-28 00:06:55 +05:30
write ( 6 , '(a24,1x,L8)' ) ' timeSyncing: ' , numerics_timeSyncing
2012-09-24 21:52:25 +05:30
write ( 6 , '(a24,1x,L8)' ) ' analytic Jacobian: ' , analyticJaco
write ( 6 , '(a24,1x,es8.1,/)' ) ' unitlength: ' , numerics_unitlength
2012-03-07 15:37:29 +05:30
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,i8)' ) ' nHomog: ' , nHomog
write ( 6 , '(a24,1x,es8.1)' ) ' subStepMinHomog: ' , subStepMinHomog
write ( 6 , '(a24,1x,es8.1)' ) ' subStepSizeHomog: ' , subStepSizeHomog
write ( 6 , '(a24,1x,es8.1)' ) ' stepIncreaseHomog: ' , stepIncreaseHomog
write ( 6 , '(a24,1x,i8,/)' ) ' nMPstate: ' , nMPstate
2009-07-31 17:32:20 +05:30
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! RGC parameters
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,es8.1)' ) ' aTol_RGC: ' , absTol_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' rTol_RGC: ' , relTol_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' aMax_RGC: ' , absMax_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' rMax_RGC: ' , relMax_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' perturbPenalty_RGC: ' , pPert_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' relevantMismatch_RGC: ' , xSmoo_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' viscosityrate_RGC: ' , viscPower_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' viscositymodulus_RGC: ' , viscModus_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' maxrelaxation_RGC: ' , maxdRelax_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' maxVolDiscrepancy_RGC: ' , maxVolDiscr_RGC
write ( 6 , '(a24,1x,es8.1)' ) ' volDiscrepancyMod_RGC: ' , volDiscrMod_RGC
write ( 6 , '(a24,1x,es8.1,/)' ) ' discrepancyPower_RGC: ' , volDiscrPow_RGC
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! Random seeding parameter
2012-06-15 21:40:21 +05:30
write ( 6 , '(a24,1x,i16,/)' ) ' fixed_seed: ' , fixedSeed
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! openMP parameter
2012-06-15 21:40:21 +05:30
!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
2010-10-13 21:34:44 +05:30
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! spectral parameters
2012-06-15 21:40:21 +05:30
#ifdef Spectral
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,es8.1)' ) ' err_div_tol: ' , err_div_tol
write ( 6 , '(a24,1x,es8.1)' ) ' err_stress_tolrel: ' , err_stress_tolrel
2012-04-11 18:27:25 +05:30
write ( 6 , '(a24,1x,es8.1)' ) ' err_stress_tolabs: ' , err_stress_tolabs
2012-08-29 00:49:47 +05:30
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,i8)' ) ' itmax: ' , itmax
write ( 6 , '(a24,1x,i8)' ) ' itmin: ' , itmin
2012-09-13 15:18:38 +05:30
write ( 6 , '(a24,1x,i8)' ) ' maxCutBack: ' , maxCutBack
2012-12-14 20:48:04 +05:30
write ( 6 , '(a24,1x,i8)' ) ' regridMode: ' , regridMode
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,L8)' ) ' memory_efficient: ' , memory_efficient
2012-03-07 15:37:29 +05:30
if ( fftw_timelimit < 0.0_pReal ) then
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,L8)' ) ' fftw_timelimit: ' , . false .
2012-03-07 15:37:29 +05:30
else
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,es8.1)' ) ' fftw_timelimit: ' , fftw_timelimit
2012-03-07 15:37:29 +05:30
endif
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,a)' ) ' fftw_plan_mode: ' , trim ( fftw_plan_mode )
2012-08-29 00:49:47 +05:30
2012-08-07 22:53:13 +05:30
write ( 6 , '(a24,1x,a)' ) ' myfilter: ' , trim ( myfilter )
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,i8)' ) ' fftw_planner_flag: ' , fftw_planner_flag
write ( 6 , '(a24,1x,es8.1)' ) ' rotation_tol: ' , rotation_tol
2013-01-16 16:10:53 +05:30
write ( 6 , '(a24,1x,i8)' ) ' divergence_correction: ' , divergence_correction
2012-03-20 23:31:31 +05:30
write ( 6 , '(a24,1x,L8,/)' ) ' update_gamma: ' , update_gamma
2012-08-29 00:49:47 +05:30
#ifdef PETSc
write ( 6 , '(a24,1x,es8.1)' ) ' err_f_tol: ' , err_f_tol
write ( 6 , '(a24,1x,es8.1)' ) ' err_p_tol: ' , err_p_tol
2013-02-28 23:05:02 +05:30
write ( 6 , '(a24,1x,es8.1)' ) ' polarAlpha: ' , polarAlpha
write ( 6 , '(a24,1x,es8.1)' ) ' polarBeta: ' , polarBeta
2012-08-29 00:49:47 +05:30
write ( 6 , '(a24,1x,a)' ) ' myspectralsolver: ' , trim ( myspectralsolver )
write ( 6 , '(a24,1x,a)' ) ' PETSc_options: ' , trim ( petsc_options )
#endif
2012-06-15 21:40:21 +05:30
#endif
2010-12-02 16:34:29 +05:30
2013-02-11 15:14:17 +05:30
!--------------------------------------------------------------------------------------------------
! sanity checks
2012-03-07 15:37:29 +05:30
if ( relevantStrain < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'relevantStrain' )
if ( defgradTolerance < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'defgradTolerance' )
if ( iJacoStiffness < 1_pInt ) call IO_error ( 301_pInt , ext_msg = 'iJacoStiffness' )
if ( iJacoLpresiduum < 1_pInt ) call IO_error ( 301_pInt , ext_msg = 'iJacoLpresiduum' )
if ( pert_Fg < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'pert_Fg' )
if ( pert_method < = 0_pInt . or . pert_method > = 4_pInt ) &
call IO_error ( 301_pInt , ext_msg = 'pert_method' )
if ( nHomog < 1_pInt ) call IO_error ( 301_pInt , ext_msg = 'nHomog' )
if ( nMPstate < 1_pInt ) call IO_error ( 301_pInt , ext_msg = 'nMPstate' )
if ( nCryst < 1_pInt ) call IO_error ( 301_pInt , ext_msg = 'nCryst' )
if ( nState < 1_pInt ) call IO_error ( 301_pInt , ext_msg = 'nState' )
if ( nStress < 1_pInt ) call IO_error ( 301_pInt , ext_msg = 'nStress' )
if ( subStepMinCryst < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'subStepMinCryst' )
if ( subStepSizeCryst < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'subStepSizeCryst' )
if ( stepIncreaseCryst < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'stepIncreaseCryst' )
if ( subStepMinHomog < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'subStepMinHomog' )
if ( subStepSizeHomog < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'subStepSizeHomog' )
if ( stepIncreaseHomog < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'stepIncreaseHomog' )
if ( rTol_crystalliteState < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'rTol_crystalliteState' )
if ( rTol_crystalliteTemperature < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'rTol_crystalliteTemperature' )
if ( rTol_crystalliteStress < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'rTol_crystalliteStress' )
if ( aTol_crystalliteStress < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'aTol_crystalliteStress' )
if ( any ( numerics_integrator < = 0_pInt ) . or . any ( numerics_integrator > = 6_pInt ) ) &
call IO_error ( 301_pInt , ext_msg = 'integrator' )
2012-09-24 21:52:25 +05:30
if ( numerics_unitlength < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'unitlength' )
2012-03-07 15:37:29 +05:30
if ( absTol_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'absTol_RGC' )
if ( relTol_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'relTol_RGC' )
if ( absMax_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'absMax_RGC' )
if ( relMax_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'relMax_RGC' )
if ( pPert_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'pPert_RGC' )
if ( xSmoo_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'xSmoo_RGC' )
if ( viscPower_RGC < 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'viscPower_RGC' )
if ( viscModus_RGC < 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'viscModus_RGC' )
if ( refRelaxRate_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'refRelaxRate_RGC' )
if ( maxdRelax_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'maxdRelax_RGC' )
if ( maxVolDiscr_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'maxVolDiscr_RGC' )
if ( volDiscrMod_RGC < 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'volDiscrMod_RGC' )
if ( volDiscrPow_RGC < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'volDiscrPw_RGC' )
2012-06-15 21:40:21 +05:30
#ifdef Spectral
2012-03-07 15:37:29 +05:30
if ( err_div_tol < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'err_div_tol' )
if ( err_stress_tolrel < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'err_stress_tolrel' )
2012-04-11 18:27:25 +05:30
if ( err_stress_tolabs < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'err_stress_tolabs' )
2012-03-07 15:37:29 +05:30
if ( itmax < = 1.0_pInt ) call IO_error ( 301_pInt , ext_msg = 'itmax' )
2012-07-24 17:51:41 +05:30
if ( itmin > itmax . or . itmin < 1_pInt ) call IO_error ( 301_pInt , ext_msg = 'itmin' )
2013-01-16 16:10:53 +05:30
if ( divergence_correction < 0_pInt . or . &
2013-02-08 21:26:58 +05:30
divergence_correction > 3_pInt ) call IO_error ( 301_pInt , ext_msg = 'divergence_correction' )
2012-09-13 15:18:38 +05:30
if ( maxCutBack < = 1.0_pInt ) call IO_error ( 301_pInt , ext_msg = 'maxCutBack' )
2012-04-11 22:58:08 +05:30
if ( update_gamma . and . &
. not . memory_efficient ) call IO_error ( error_ID = 847_pInt )
2012-08-29 00:49:47 +05:30
#ifdef PETSc
if ( err_f_tol < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'err_f_tol' )
if ( err_p_tol < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'err_p_tol' )
2013-02-28 23:05:02 +05:30
if ( polarAlpha < = 0.0_pReal . or . &
polarAlpha > 2.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'polarAlpha' )
if ( polarBeta < = 0.0_pReal . or . &
polarBeta > 2.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'polarBeta' )
if ( err_p_tol < = 0.0_pReal ) call IO_error ( 301_pInt , ext_msg = 'err_p_tol' )
2012-08-29 00:49:47 +05:30
#endif
2012-06-15 21:40:21 +05:30
#endif
2013-02-11 15:14:17 +05:30
if ( fixedSeed < = 0_pInt ) &
write ( 6 , '(a,/)' ) ' No fixed Seed: Random is random!'
2012-03-07 15:37:29 +05:30
end subroutine numerics_init
2009-06-15 18:41:21 +05:30
2012-03-07 15:37:29 +05:30
end module numerics