diff --git a/VERSION b/VERSION index db8be18b4..e64a6d9e6 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0-alpha7-558-gf6c7fd280 +3.0.0-alpha7-584-gaf7452edd diff --git a/src/CLI.f90 b/src/CLI.f90 index ac1353b10..011b80d59 100644 --- a/src/CLI.f90 +++ b/src/CLI.f90 @@ -16,6 +16,7 @@ module CLI use prec use parallelization use system_routines + use IO implicit none(type,external) private @@ -23,7 +24,8 @@ module CLI CLI_restartInc = 0 !< Increment at which calculation starts character(len=:), allocatable, public, protected :: & CLI_geomFile, & !< parameter given for geometry file - CLI_loadFile !< parameter given for load case file + CLI_loadFile, & !< parameter given for load case file + CLI_materialFile public :: & getSolverJobName, & @@ -35,158 +37,187 @@ contains !> @brief initializes the solver by interpreting the command line arguments. Also writes !! information on computation to screen !-------------------------------------------------------------------------------------------------- -subroutine CLI_init +subroutine CLI_init() #include #if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINORPETSC_MINOR_MAX -- UNSUPPORTED PETSc VERSION --- UNSUPPORTED PETSc VERSION --- UNSUPPORTED PETSc VERSION --- #endif - character(len=pPathLen*3+pSTRLEN) :: & - commandLine !< command line call as string - character(len=pPathLen) :: & + character(len=:), allocatable :: & + commandLine, & !< command line call as string arg, & !< individual argument - loadCaseArg = '', & !< -l argument given to the executable - geometryArg = '', & !< -g argument given to the executable - workingDirArg = '' !< -w argument given to the executable + loadCaseArg, & !< -l argument given to the executable + geometryArg, & !< -g argument given to the executable + materialArg, & !< -m argument given to the executable + workingDirArg !< -w argument given to the executable integer :: & stat, & i integer, dimension(8) :: & dateAndTime - integer :: err external :: & quit + workingDirArg = getCWD() + print'(/,1x,a)', '<<<+- CLI init -+>>>' ! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK%203 #ifdef DEBUG print*, achar(27)//'[31m' - print'(a,/)', ' debug version - debug version - debug version - debug version - debug version' + print'(1x,a,/)', 'debug version - debug version - debug version - debug version - debug version' #else - print*, achar(27)//'[94m' + print '(a)', achar(27)//'[94m' #endif - print*, ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/ _/_/_/' - print*, ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/ _/' - print*, ' _/ _/ _/_/_/_/ _/ _/ _/ _/_/_/_/ _/_/ _/_/ _/_/' - print*, ' _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/' - print*, ' _/_/_/ _/ _/ _/ _/ _/ _/ _/_/_/ _/ _/ _/_/_/' + print '(1x,a)', ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/ _/_/_/' + print '(1x,a)', ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/ _/' + print '(1x,a)', ' _/ _/ _/_/_/_/ _/ _/ _/ _/_/_/_/ _/_/ _/_/ _/_/' + print '(1x,a)', ' _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/' + print '(1x,a)', '_/_/_/ _/ _/ _/ _/ _/ _/ _/_/_/ _/ _/ _/_/_/' #if defined(GRID) - print*, ' Grid solver' + print '(1x,a)', 'Grid solver' #elif defined(MESH) - print*, ' Mesh solver' + print '(1x,a)', 'Mesh solver' #endif #ifdef DEBUG - print'(/,a)', ' debug version - debug version - debug version - debug version - debug version' + print'(/,1x,a)', 'debug version - debug version - debug version - debug version - debug version' #endif - print*, achar(27)//'[0m' + print '(a)', achar(27)//'[0m' - print*, 'F. Roters et al., Computational Materials Science 158:420–478, 2019' - print*, 'https://doi.org/10.1016/j.commatsci.2018.04.030' + print '(1x,a)', 'F. Roters et al., Computational Materials Science 158:420–478, 2019' + print '(1x,a)', 'https://doi.org/10.1016/j.commatsci.2018.04.030' - print'(/,a)', ' Version: '//DAMASKVERSION + print '(/,1x,a)', 'Version: '//DAMASKVERSION - print'(/,a)', ' Compiled with: '//compiler_version() - print'(a)', ' Compiled on: '//CMAKE_SYSTEM - print'(a)', ' Compiler options: '//compiler_options() + print '(/,1x,a)', 'Compiled with: '//compiler_version() + print '(1x,a)', 'Compiled on: '//CMAKE_SYSTEM + print '(1x,a)', 'Compiler options: '//compiler_options() ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md - print'(/,a)', ' Compiled on: '//__DATE__//' at '//__TIME__ + print '(/,1x,a)', 'Compiled on: '//__DATE__//' at '//__TIME__ - print'(/,a,i0,a,i0,a,i0)', & - ' PETSc version: ',PETSC_VERSION_MAJOR,'.',PETSC_VERSION_MINOR,'.',PETSC_VERSION_SUBMINOR + print '(/,1x,a,1x,i0,a,i0,a,i0)', & + 'PETSc version:',PETSC_VERSION_MAJOR,'.',PETSC_VERSION_MINOR,'.',PETSC_VERSION_SUBMINOR call date_and_time(values = dateAndTime) - print'(/,a,2(i2.2,a),i4.4)', ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1) - print'(a,2(i2.2,a),i2.2)', ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7) + print '(/,1x,a,1x,2(i2.2,a),i4.4)', 'Date:',dateAndTime(3),'/',dateAndTime(2),'/',dateAndTime(1) + print '(1x,a,1x,2(i2.2,a),i2.2)', 'Time:',dateAndTime(5),':',dateAndTime(6),':',dateAndTime(7) do i = 1, command_argument_count() - call get_command_argument(i,arg,status=err) - if (err /= 0) call quit(1) + arg = getArg(i) select case(trim(arg)) ! extract key case ('-h','--help') - print'(/,a)',' #######################################################################' - print'(a)', ' DAMASK Command Line Interface:' - print'(a)', ' Düsseldorf Advanced Material Simulation Kit with PETSc-based solvers' - print'(a,/)',' #######################################################################' - print'(a,/)',' Valid command line switches:' - print'(a)', ' --geom (-g, --geometry)' - print'(a)', ' --load (-l, --loadcase)' - print'(a)', ' --workingdir (-w, --wd, --workingdirectory)' - print'(a)', ' --restart (-r, --rs)' - print'(a)', ' --help (-h)' - print'(/,a)',' -----------------------------------------------------------------------' - print'(a)', ' Mandatory arguments:' - print'(/,a)',' --geom PathToGeomFile/NameOfGeom' - print'(a)', ' Specifies the location of the geometry definition file.' - print'(/,a)',' --load PathToLoadFile/NameOfLoadFile' - print'(a)', ' Specifies the location of the load case definition file.' - print'(/,a)',' -----------------------------------------------------------------------' - print'(a)', ' Optional arguments:' - print'(/,a)',' --workingdirectory PathToWorkingDirectory' - print'(a)', ' Specifies the working directory and overwrites the default ./' - print'(a)', ' Make sure the file "material.yaml" exists in the working' - print'(a)', ' directory.' - print'(a)', ' For further configuration place "numerics.yaml"' - print'(a)',' in that directory.' - print'(/,a)',' --restart N' - print'(a)', ' Reads in increment N and continues with calculating' - print'(a)', ' increment N+1, N+2, ... based on this.' - print'(a)', ' Appends to existing results file' - print'(a)', ' "NameOfGeom_NameOfLoadFile.hdf5".' - print'(a)', ' Works only if the restart information for increment N' - print'(a)', ' is available in the working directory.' - print'(/,a)',' -----------------------------------------------------------------------' - print'(a)', ' Help:' - print'(/,a)',' --help' - print'(a,/)',' Prints this message and exits' + print '(/,1x,a)','#######################################################################' + print '(1x,a)', 'DAMASK Command Line Interface:' + print '(1x,a)', 'Düsseldorf Advanced Material Simulation Kit with PETSc-based solvers' + print '(1x,a,/)','#######################################################################' + print '(1x,a,/)','Valid command line switches:' + print '(1x,a)', ' --geom (-g, --geometry)' + print '(1x,a)', ' --load (-l, --loadcase)' + print '(1x,a)', ' --material (-m, --materialconfig)' + print '(1x,a)', ' --workingdir (-w, --wd, --workingdirectory)' + print '(1x,a)', ' --restart (-r, --rs)' + print '(1x,a)', ' --help (-h)' + print '(/,1x,a)','-----------------------------------------------------------------------' + print '(1x,a)', 'Mandatory arguments:' + print '(/,1x,a)',' --geom PathToGeomFile/NameOfGeom' + print '(1x,a)', ' Specifies the location of the geometry definition file.' + print '(/,1x,a)',' --load PathToLoadFile/NameOfLoadFile' + print '(1x,a)', ' Specifies the location of the load case definition file.' + print '(/,1x,a)',' --material PathToMaterialConfigurationFile/NameOfMaterialConfigurationFile' + print '(1x,a)', ' Specifies the location of the material configuration file.' + print '(/,1x,a)','-----------------------------------------------------------------------' + print '(1x,a)', 'Optional arguments:' + print '(/,1x,a)',' --workingdirectory PathToWorkingDirectory' + print '(1x,a)', ' Specifies the base directory of relative paths.' + print '(/,1x,a)',' --restart N' + print '(1x,a)', ' Reads in increment N and continues with calculating' + print '(1x,a)', ' increment N+1, N+2, ... based on this.' + print '(1x,a)', ' Appends to existing results file' + print '(1x,a)', ' "NameOfGeom_NameOfLoadFile_NameOfMaterialConfigurationFile.hdf5".' + print '(1x,a)', ' Works only if the restart information for increment N' + print '(1x,a)', ' is available in the base directory.' + print '(/,1x,a)','-----------------------------------------------------------------------' + print '(1x,a)', 'Help:' + print '(/,1x,a)',' --help' + print '(1x,a,/)',' Prints this message and exits' call quit(0) ! normal Termination case ('-l', '--load', '--loadcase') - call get_command_argument(i+1,loadCaseArg,status=err) + loadCaseArg = getArg(i+1) case ('-g', '--geom', '--geometry') - call get_command_argument(i+1,geometryArg,status=err) + geometryArg = getArg(i+1) + case ('-m', '--material', '--materialconfig') + materialArg = getArg(i+1) case ('-w', '--wd', '--workingdir', '--workingdirectory') - call get_command_argument(i+1,workingDirArg,status=err) + workingDirArg = getArg(i+1) case ('-r', '--rs', '--restart') - call get_command_argument(i+1,arg,status=err) + arg = getArg(i+1) read(arg,*,iostat=stat) CLI_restartInc - if (CLI_restartInc < 0 .or. stat /=0) then - print'(/,a)', ' ERROR: Could not parse restart increment: '//trim(arg) + if (CLI_restartInc < 0 .or. stat /= 0) then + print'(/,1x,a)', 'ERROR: Could not parse restart increment: '//trim(arg) call quit(1) end if end select - if (err /= 0) call quit(1) end do - if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then - print'(/,a)', ' ERROR: Please specify geometry AND load case (-h for help)' + if (.not. all([allocated(loadcaseArg),allocated(geometryArg),allocated(materialArg)])) then + print'(/,1x,a)', 'ERROR: Please specify geometry AND load case AND material configuration (-h for help)' call quit(1) end if - if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg)) - CLI_geomFile = getGeometryFile(geometryArg) - CLI_loadFile = getLoadCaseFile(loadCaseArg) + call setWorkingDirectory(trim(workingDirArg)) + CLI_geomFile = getPathRelCWD(geometryArg,'geometry') + CLI_loadFile = getPathRelCWD(loadCaseArg,'load case') + CLI_materialFile = getPathRelCWD(materialArg,'material configuration') - call get_command(commandLine) - print'(/,a)', ' Host name: '//getHostName() - print'(a)', ' User name: '//getUserName() + commandLine = getArg(-1) - print'(/a)', ' Command line call: '//trim(commandLine) - if (len_trim(workingDirArg) > 0) & - print'(a)', ' Working dir argument: '//trim(workingDirArg) - print'(a)', ' Geometry argument: '//trim(geometryArg) - print'(a)', ' Load case argument: '//trim(loadcaseArg) - print'(/,a)', ' Working directory: '//getCWD() - print'(a)', ' Geometry file: '//CLI_geomFile - print'(a)', ' Load case file: '//CLI_loadFile - print'(a)', ' Solver job name: '//getSolverJobName() + print'(/,1x,a)', 'Host name: '//getHostName() + print'(1x,a)', 'User name: '//getUserName() + + print'(/,1x,a,/)', 'Command line call: '//trim(commandLine) + print'(1x,a)', 'Working directory: '//IO_glueDiffering(getCWD(),workingDirArg) + print'(1x,a)', 'Geometry: '//IO_glueDiffering(CLI_geomFile,geometryArg) + print'(1x,a)', 'Load case: '//IO_glueDiffering(CLI_loadFile,loadCaseArg) + print'(1x,a)', 'Material config: '//IO_glueDiffering(CLI_materialFile,materialArg) + print'(1x,a)', 'Solver job name: '//getSolverJobName() if (CLI_restartInc > 0) & - print'(a,i6.6)', ' Restart from increment: ', CLI_restartInc + print'(1x,a,i6.6)', 'Restart from increment: ', CLI_restartInc + end subroutine CLI_init +!-------------------------------------------------------------------------------------------------- +!> @brief Get argument from command line. +!-------------------------------------------------------------------------------------------------- +function getArg(n) + + integer, intent(in) :: n !< number of the argument + character(len=:), allocatable :: getArg + + integer :: l,err + external :: quit + + + allocate(character(len=0)::getArg) + if (n<0) then + call get_command(getArg, length=l) + else + call get_command_argument(n,getArg,length=l) + endif + deallocate(getArg) + allocate(character(len=l)::getArg) + if (n<0) then + call get_command(getArg, status=err) + else + call get_command_argument(n,getArg,status=err) + endif + if (err /= 0) call quit(1) + +end function getArg + !-------------------------------------------------------------------------------------------------- !> @brief extract working directory from given argument or from location of geometry file, @@ -196,9 +227,11 @@ subroutine setWorkingDirectory(workingDirectoryArg) character(len=*), intent(in) :: workingDirectoryArg !< working directory argument character(len=:), allocatable :: workingDirectory + logical :: error external :: quit + absolutePath: if (workingDirectoryArg(1:1) == '/') then workingDirectory = workingDirectoryArg else absolutePath @@ -206,10 +239,10 @@ subroutine setWorkingDirectory(workingDirectoryArg) workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg end if absolutePath - workingDirectory = trim(rectifyPath(workingDirectory)) + workingDirectory = trim(normpath(workingDirectory)) error = setCWD(trim(workingDirectory)) if (error) then - print*, 'ERROR: Invalid Working directory: '//trim(workingDirectory) + print '(1x,a)', 'ERROR: Invalid Working directory: '//trim(workingDirectory) call quit(1) end if @@ -222,8 +255,10 @@ end subroutine setWorkingDirectory function getSolverJobName() character(len=:), allocatable :: getSolverJobName + integer :: posExt,posSep + posExt = scan(CLI_geomFile,'.',back=.true.) posSep = scan(CLI_geomFile,'/',back=.true.) @@ -238,123 +273,107 @@ end function getSolverJobName !-------------------------------------------------------------------------------------------------- -!> @brief basename of geometry file with extension from command line arguments +!> @brief Translate path as relative to CWD and check for existence. !-------------------------------------------------------------------------------------------------- -function getGeometryFile(geometryParameter) +function getPathRelCWD(path,fileType) + + character(len=:), allocatable :: getPathRelCWD + character(len=*), intent(in) :: path + character(len=*), intent(in) :: fileType - character(len=:), allocatable :: getGeometryFile - character(len=*), intent(in) :: geometryParameter logical :: file_exists external :: quit - getGeometryFile = trim(geometryParameter) - if (scan(getGeometryFile,'/') /= 1) getGeometryFile = getCWD()//'/'//trim(getGeometryFile) - getGeometryFile = trim(makeRelativePath(getCWD(), getGeometryFile)) - inquire(file=getGeometryFile, exist=file_exists) + getPathRelCWD = trim(path) + if (scan(getPathRelCWD,'/') /= 1) getPathRelCWD = getCWD()//'/'//trim(getPathRelCWD) + getPathRelCWD = trim(relpath(getPathRelCWD,getCWD())) + + inquire(file=getPathRelCWD, exist=file_exists) if (.not. file_exists) then - print*, 'ERROR: Geometry file does not exists: '//trim(getGeometryFile) + print '(1x,a)', 'ERROR: '//fileType//' file does not exist: '//trim(getPathRelCWD) call quit(1) end if -end function getGeometryFile +end function getPathRelCWD !-------------------------------------------------------------------------------------------------- -!> @brief relative path of load case from command line arguments +!> @brief Remove ../, /./, and // from path. +!> @details Works only if absolute path is given. !-------------------------------------------------------------------------------------------------- -function getLoadCaseFile(loadCaseParameter) - - character(len=:), allocatable :: getLoadCaseFile - character(len=*), intent(in) :: loadCaseParameter - logical :: file_exists - external :: quit - - getLoadCaseFile = trim(loadCaseParameter) - if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = getCWD()//'/'//trim(getLoadCaseFile) - getLoadCaseFile = trim(makeRelativePath(getCWD(), getLoadCaseFile)) - - inquire(file=getLoadCaseFile, exist=file_exists) - if (.not. file_exists) then - print*, 'ERROR: Load case file does not exists: '//trim(getLoadCaseFile) - call quit(1) - end if - -end function getLoadCaseFile - - -!-------------------------------------------------------------------------------------------------- -!> @brief remove ../, /./, and // from path. -!> @details works only if absolute path is given -!-------------------------------------------------------------------------------------------------- -function rectifyPath(path) +function normpath(path) character(len=*), intent(in) :: path - character(len=:), allocatable :: rectifyPath + character(len=:), allocatable :: normpath + integer :: i,j,k,l + !-------------------------------------------------------------------------------------------------- ! remove /./ from path - rectifyPath = trim(path) - l = len_trim(rectifyPath) + normpath = trim(path) + l = len_trim(normpath) do i = l,3,-1 - if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' + if (normpath(i-2:i) == '/./') normpath(i-1:l) = normpath(i+1:l)//' ' end do !-------------------------------------------------------------------------------------------------- ! remove // from path - l = len_trim(rectifyPath) + l = len_trim(normpath) do i = l,2,-1 - if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' ' + if (normpath(i-1:i) == '//') normpath(i-1:l) = normpath(i:l)//' ' end do !-------------------------------------------------------------------------------------------------- -! remove ../ and corresponding directory from rectifyPath - l = len_trim(rectifyPath) - i = index(rectifyPath(i:l),'../') +! remove ../ and corresponding directory from path + l = len_trim(normpath) + i = index(normpath(i:l),'../') j = 0 do while (i > j) - j = scan(rectifyPath(1:i-2),'/',back=.true.) - rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j) - if (rectifyPath(j+1:j+1) == '/') then !search for '//' that appear in case of XXX/../../XXX - k = len_trim(rectifyPath) - rectifyPath(j+1:k-1) = rectifyPath(j+2:k) - rectifyPath(k:k) = ' ' + j = scan(normpath(1:i-2),'/',back=.true.) + normpath(j+1:l) = normpath(i+3:l)//repeat(' ',2+i-j) + if (normpath(j+1:j+1) == '/') then !search for '//' that appear in case of XXX/../../XXX + k = len_trim(normpath) + normpath(j+1:k-1) = normpath(j+2:k) + normpath(k:k) = ' ' end if - i = j+index(rectifyPath(j+1:l),'../') + i = j+index(normpath(j+1:l),'../') end do - if (len_trim(rectifyPath) == 0) rectifyPath = '/' + if (len_trim(normpath) == 0) normpath = '/' - rectifyPath = trim(rectifyPath) + normpath = trim(normpath) -end function rectifyPath +end function normpath !-------------------------------------------------------------------------------------------------- -!> @brief Determine relative path from absolute a to absolute b +!> @brief Determine relative path. !-------------------------------------------------------------------------------------------------- -function makeRelativePath(a,b) +function relpath(path,start) - character(len=*), intent(in) :: a,b - character(len=pPathLen) :: a_cleaned,b_cleaned - character(len=:), allocatable :: makeRelativePath + character(len=*), intent(in) :: start,path + character(len=:), allocatable :: relpath + + character(len=:), allocatable :: start_cleaned,path_cleaned integer :: i,posLastCommonSlash,remainingSlashes + posLastCommonSlash = 0 remainingSlashes = 0 - a_cleaned = rectifyPath(trim(a)//'/') - b_cleaned = rectifyPath(b) + start_cleaned = normpath(trim(start)//'/') + path_cleaned = normpath(path) - do i = 1, min(len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned))) - if (a_cleaned(i:i) /= b_cleaned(i:i)) exit - if (a_cleaned(i:i) == '/') posLastCommonSlash = i + do i = 1, min(len_trim(start_cleaned),len_trim(path_cleaned)) + if (start_cleaned(i:i) /= path_cleaned(i:i)) exit + if (start_cleaned(i:i) == '/') posLastCommonSlash = i end do - do i = posLastCommonSlash+1,len_trim(a_cleaned) - if (a_cleaned(i:i) == '/') remainingSlashes = remainingSlashes + 1 + do i = posLastCommonSlash+1,len_trim(start_cleaned) + if (start_cleaned(i:i) == '/') remainingSlashes = remainingSlashes + 1 end do - makeRelativePath = repeat('..'//'/',remainingSlashes)//b_cleaned(posLastCommonSlash+1:len_trim(b_cleaned)) + relpath = repeat('..'//'/',remainingSlashes)//path_cleaned(posLastCommonSlash+1:len_trim(path_cleaned)) -end function makeRelativePath +end function relpath end module CLI diff --git a/src/IO.f90 b/src/IO.f90 index 27e650825..39a48e1e5 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -38,6 +38,7 @@ module IO IO_realValue, & IO_lc, & IO_rmComment, & + IO_glueDiffering, & IO_intAsStr, & IO_strAsInt, & IO_strAsReal, & @@ -319,6 +320,7 @@ function IO_rmComment(line) character(len=*), intent(in) :: line character(len=:), allocatable :: IO_rmComment + integer :: split @@ -333,15 +335,35 @@ function IO_rmComment(line) end function IO_rmComment +!-------------------------------------------------------------------------------------------------- +! @brief Return first (with glued on second if they differ). +!-------------------------------------------------------------------------------------------------- +function IO_glueDiffering(first,second,glue) + + character(len=*), intent(in) :: first + character(len=*), intent(in) :: second + character(len=*), optional, intent(in) :: glue + character(len=:), allocatable :: IO_glueDiffering + + character(len=:), allocatable :: glue_ + + + glue_ = misc_optional(glue,'<--') + IO_glueDiffering = trim(first) + if (trim(first) /= trim(second)) IO_glueDiffering = IO_glueDiffering//' '//trim(glue_)//' '//trim(second) + +end function IO_glueDiffering + + !-------------------------------------------------------------------------------------------------- !> @brief Return given int value as string. !-------------------------------------------------------------------------------------------------- function IO_intAsStr(i) integer, intent(in) :: i - character(len=:), allocatable :: IO_intAsStr + allocate(character(len=merge(2,1,i<0) + floor(log10(real(abs(merge(1,i,i==0))))))::IO_intAsStr) write(IO_intAsStr,'(i0)') i diff --git a/src/config.f90 b/src/config.f90 index 6e173db57..9006839bf 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -9,7 +9,9 @@ module config use YAML_types use result use parallelization - +#if defined(MESH) || defined(GRID) + use CLI +#endif implicit none(type,external) private @@ -96,17 +98,20 @@ end function config_listReferences subroutine parse_material() logical :: fileExists - character(len=:), allocatable :: fileContent - - - inquire(file='material.yaml',exist=fileExists) - if (.not. fileExists) call IO_error(100,ext_msg='material.yaml') + character(len=:), allocatable :: & + fileContent, fname if (worldrank == 0) then - print'(/,1x,a)', 'reading material.yaml'; flush(IO_STDOUT) - fileContent = IO_read('material.yaml') + print'(/,1x,a)', 'reading material configuration'; flush(IO_STDOUT) +#if defined(MESH) || defined(GRID) + fname = CLI_materialFile +#else + fname = 'material.yaml' +#endif + fileContent = IO_read(fname) + if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) call result_openJobFile(parallel=.false.) - call result_writeDataset_str(fileContent,'setup','material.yaml','main configuration') + call result_writeDataset_str(fileContent,'setup',fname,'material configuration') call result_closeJobFile() end if call parallelization_bcast_str(fileContent) diff --git a/src/lattice.f90 b/src/lattice.f90 index 350860ecb..c44e775db 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -40,19 +40,19 @@ module lattice real(pREAL), dimension(3+3,CF_NSLIP), parameter :: & CF_SYSTEMSLIP = reshape(real([& - ! <110>{111} systems - 0, 1,-1, 1, 1, 1, & ! B2 - -1, 0, 1, 1, 1, 1, & ! B4 - 1,-1, 0, 1, 1, 1, & ! B5 - 0,-1,-1, -1,-1, 1, & ! C1 - 1, 0, 1, -1,-1, 1, & ! C3 - -1, 1, 0, -1,-1, 1, & ! C5 - 0,-1, 1, 1,-1,-1, & ! A2 - -1, 0,-1, 1,-1,-1, & ! A3 - 1, 1, 0, 1,-1,-1, & ! A6 - 0, 1, 1, -1, 1,-1, & ! D1 - 1, 0,-1, -1, 1,-1, & ! D4 - -1,-1, 0, -1, 1,-1, & ! D6 + ! <110>{111} systems (Thompson tetrahedron labeling according to Fig. 3 of 10.1016/S1572-4859(05)80003-8) + 0, 1,-1, 1, 1, 1, & ! AC(d) + -1, 0, 1, 1, 1, 1, & ! CB(d) + 1,-1, 0, 1, 1, 1, & ! BA(d) + 0,-1,-1, -1,-1, 1, & ! BD(c) + 1, 0, 1, -1,-1, 1, & ! DA(c) + -1, 1, 0, -1,-1, 1, & ! AB(c) + 0,-1, 1, 1,-1,-1, & ! CA(b) + -1, 0,-1, 1,-1,-1, & ! AD(b) + 1, 1, 0, 1,-1,-1, & ! DC(b) + 0, 1, 1, -1, 1,-1, & ! DB(a) + 1, 0,-1, -1, 1,-1, & ! BC(a) + -1,-1, 0, -1, 1,-1, & ! CD(a) ! <110>{110}/non-octahedral systems 1, 1, 0, 1,-1, 0, & 1,-1, 0, 1, 1, 0, & @@ -123,31 +123,31 @@ module lattice real(pREAL), dimension(3+3,CI_NSLIP), parameter :: & CI_SYSTEMSLIP = reshape(real([& ! <111>{110} systems - 1,-1, 1, 0, 1, 1, & ! D1 - -1,-1, 1, 0, 1, 1, & ! C1 - 1, 1, 1, 0,-1, 1, & ! B2 - -1, 1, 1, 0,-1, 1, & ! A2 - -1, 1, 1, 1, 0, 1, & ! A3 - -1,-1, 1, 1, 0, 1, & ! C3 - 1, 1, 1, -1, 0, 1, & ! B4 - 1,-1, 1, -1, 0, 1, & ! D4 - -1, 1, 1, 1, 1, 0, & ! A6 - -1, 1,-1, 1, 1, 0, & ! D6 - 1, 1, 1, -1, 1, 0, & ! B5 - 1, 1,-1, -1, 1, 0, & ! C5 + 1,-1, 1, 0, 1, 1, & + -1,-1, 1, 0, 1, 1, & + 1, 1, 1, 0,-1, 1, & + -1, 1, 1, 0,-1, 1, & + -1, 1, 1, 1, 0, 1, & + -1,-1, 1, 1, 0, 1, & + 1, 1, 1, -1, 0, 1, & + 1,-1, 1, -1, 0, 1, & + -1, 1, 1, 1, 1, 0, & + -1, 1,-1, 1, 1, 0, & + 1, 1, 1, -1, 1, 0, & + 1, 1,-1, -1, 1, 0, & ! <111>{112} systems - -1, 1, 1, 2, 1, 1, & ! A-4 - 1, 1, 1, -2, 1, 1, & ! B-3 - 1, 1,-1, 2,-1, 1, & ! C-10 - 1,-1, 1, 2, 1,-1, & ! D-9 - 1,-1, 1, 1, 2, 1, & ! D-6 - 1, 1,-1, -1, 2, 1, & ! C-5 - 1, 1, 1, 1,-2, 1, & ! B-12 - -1, 1, 1, 1, 2,-1, & ! A-11 - 1, 1,-1, 1, 1, 2, & ! C-2 - 1,-1, 1, -1, 1, 2, & ! D-1 - -1, 1, 1, 1,-1, 2, & ! A-8 - 1, 1, 1, 1, 1,-2, & ! B-7 + -1, 1, 1, 2, 1, 1, & + 1, 1, 1, -2, 1, 1, & + 1, 1,-1, 2,-1, 1, & + 1,-1, 1, 2, 1,-1, & + 1,-1, 1, 1, 2, 1, & + 1, 1,-1, -1, 2, 1, & + 1, 1, 1, 1,-2, 1, & + -1, 1, 1, 1, 2,-1, & + 1, 1,-1, 1, 1, 2, & + 1,-1, 1, -1, 1, 2, & + -1, 1, 1, 1,-1, 2, & + 1, 1, 1, 1, 1,-2, & ! Slip system <111>{123} 1, 1,-1, 1, 2, 3, & 1,-1, 1, -1, 2, 3, & diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index eff65f9f3..be859dddb 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -82,13 +82,14 @@ module function plastic_isotropic_init() result(myPlasticity) do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), stt => state(ph)) + associate(prm => param(ph), & + stt => state(ph)) phase => phases%get_dict(ph) mech => phase%get_dict('mechanical') pl => mech%get_dict('plastic') - print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph) + print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph) refs = config_listReferences(pl,indent=3) if (len(refs) > 0) print'(/,1x,a)', refs diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index ad2543b83..0bebe3368 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -116,14 +116,15 @@ module function plastic_kinehardening_init() result(myPlasticity) do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph), & + associate(prm => param(ph), & + stt => state(ph), dlt => deltaState(ph), & idx_dot => indexDotState(ph)) phase => phases%get_dict(ph) mech => phase%get_dict('mechanical') pl => mech%get_dict('plastic') - print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph) + print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph) refs = config_listReferences(pl,indent=3) if (len(refs) > 0) print'(/,1x,a)', refs @@ -150,28 +151,21 @@ module function plastic_kinehardening_init() result(myPlasticity) prm%P_nS_pos = prm%P prm%P_nS_neg = prm%P end if - prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), & - phase_lattice(ph)) - - xi_0 = pl%get_as1dReal('xi_0', requiredSize=size(N_sl)) - prm%xi_inf = pl%get_as1dReal('xi_inf', requiredSize=size(N_sl)) - prm%chi_inf = pl%get_as1dReal('chi_inf', requiredSize=size(N_sl)) - prm%h_0_xi = pl%get_as1dReal('h_0_xi', requiredSize=size(N_sl)) - prm%h_0_chi = pl%get_as1dReal('h_0_chi', requiredSize=size(N_sl)) - prm%h_inf_xi = pl%get_as1dReal('h_inf_xi', requiredSize=size(N_sl)) - prm%h_inf_chi = pl%get_as1dReal('h_inf_chi', requiredSize=size(N_sl)) prm%dot_gamma_0 = pl%get_asReal('dot_gamma_0') prm%n = pl%get_asReal('n') - ! expand: family => system - xi_0 = math_expand(xi_0, N_sl) - prm%xi_inf = math_expand(prm%xi_inf, N_sl) - prm%chi_inf = math_expand(prm%chi_inf, N_sl) - prm%h_0_xi = math_expand(prm%h_0_xi, N_sl) - prm%h_0_chi = math_expand(prm%h_0_chi, N_sl) - prm%h_inf_xi = math_expand(prm%h_inf_xi, N_sl) - prm%h_inf_chi = math_expand(prm%h_inf_chi, N_sl) + prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), & + phase_lattice(ph)) + + xi_0 = math_expand(pl%get_as1dReal('xi_0', requiredSize=size(N_sl)),N_sl) + prm%xi_inf = math_expand(pl%get_as1dReal('xi_inf', requiredSize=size(N_sl)),N_sl) + prm%chi_inf = math_expand(pl%get_as1dReal('chi_inf', requiredSize=size(N_sl)),N_sl) + prm%h_0_xi = math_expand(pl%get_as1dReal('h_0_xi', requiredSize=size(N_sl)),N_sl) + prm%h_0_chi = math_expand(pl%get_as1dReal('h_0_chi', requiredSize=size(N_sl)),N_sl) + prm%h_inf_xi = math_expand(pl%get_as1dReal('h_inf_xi', requiredSize=size(N_sl)),N_sl) + prm%h_inf_chi = math_expand(pl%get_as1dReal('h_inf_chi', requiredSize=size(N_sl)),N_sl) + !-------------------------------------------------------------------------------------------------- ! sanity checks @@ -183,7 +177,13 @@ module function plastic_kinehardening_init() result(myPlasticity) else slipActive xi_0 = emptyRealArray - allocate(prm%xi_inf,prm%chi_inf,prm%h_0_xi,prm%h_inf_xi,prm%h_0_chi,prm%h_inf_chi,source=emptyRealArray) + allocate(prm%xi_inf, & + prm%chi_inf, & + prm%h_0_xi, & + prm%h_0_chi, & + prm%h_inf_xi, & + prm%h_inf_chi, & + source=emptyRealArray) allocate(prm%h_sl_sl(0,0)) end if slipActive diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 index b79a61183..748e0579b 100644 --- a/src/phase_mechanical_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -29,9 +29,12 @@ module function plastic_none_init() result(myPlasticity) phases => config_material%get_dict('phase') + do ph = 1, phases%length if (.not. myPlasticity(ph)) cycle - print'(a,i0,a)', ' phase ',ph + + print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph) + call phase_allocateState(plasticState(ph),count(material_ID_phase == ph),0,0,0) end do diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index d50f562ca..4e45066b5 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -92,17 +92,13 @@ submodule(phase:plastic) nonlocal mu, & nu real(pREAL), dimension(:), allocatable :: & - d_ed, & !< minimum stable edge dipole height - d_sc, & !< minimum stable screw dipole height - tau_Peierls_ed, & - tau_Peierls_sc, & i_sl, & !< mean free path prefactor for each b_sl !< absolute length of Burgers vector [m] real(pREAL), dimension(:,:), allocatable :: & slip_normal, & slip_direction, & slip_transverse, & - minDipoleHeight, & ! edge and screw + minDipoleHeight, & ! minimum stable dipole height edge and screw peierlsstress, & ! edge and screw h_sl_sl ,& !< coefficients for slip-slip interaction forestProjection_Edge, & !< matrix of forest projections of edge dislocations @@ -236,7 +232,7 @@ module function plastic_nonlocal_init() result(myPlasticity) mech => phase%get_dict('mechanical') pl => mech%get_dict('plastic') - print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph) + print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph) refs = config_listReferences(pl,indent=3) if (len(refs) > 0) print'(/,1x,a)', refs @@ -295,52 +291,41 @@ module function plastic_nonlocal_init() result(myPlasticity) ini%rho_d_ed_0 = pl%get_as1dReal('rho_d_ed_0', requiredSize=size(ini%N_sl)) ini%rho_d_sc_0 = pl%get_as1dReal('rho_d_sc_0', requiredSize=size(ini%N_sl)) - prm%i_sl = pl%get_as1dReal('i_sl', requiredSize=size(ini%N_sl)) - prm%b_sl = pl%get_as1dReal('b_sl', requiredSize=size(ini%N_sl)) + prm%i_sl = math_expand(pl%get_as1dReal('i_sl', requiredSize=size(ini%N_sl)),ini%N_sl) + prm%b_sl = math_expand(pl%get_as1dReal('b_sl', requiredSize=size(ini%N_sl)),ini%N_sl) - prm%i_sl = math_expand(prm%i_sl,ini%N_sl) - prm%b_sl = math_expand(prm%b_sl,ini%N_sl) - - prm%d_ed = pl%get_as1dReal('d_ed', requiredSize=size(ini%N_sl)) - prm%d_sc = pl%get_as1dReal('d_sc', requiredSize=size(ini%N_sl)) - prm%d_ed = math_expand(prm%d_ed,ini%N_sl) - prm%d_sc = math_expand(prm%d_sc,ini%N_sl) allocate(prm%minDipoleHeight(prm%sum_N_sl,2)) - prm%minDipoleHeight(:,1) = prm%d_ed - prm%minDipoleHeight(:,2) = prm%d_sc + prm%minDipoleHeight(:,1) = math_expand(pl%get_as1dReal('d_ed', requiredSize=size(ini%N_sl)),ini%N_sl) + prm%minDipoleHeight(:,2) = math_expand(pl%get_as1dReal('d_sc', requiredSize=size(ini%N_sl)),ini%N_sl) - prm%tau_Peierls_ed = pl%get_as1dReal('tau_Peierls_ed', requiredSize=size(ini%N_sl)) - prm%tau_Peierls_sc = pl%get_as1dReal('tau_Peierls_sc', requiredSize=size(ini%N_sl)) - prm%tau_Peierls_ed = math_expand(prm%tau_Peierls_ed,ini%N_sl) - prm%tau_Peierls_sc = math_expand(prm%tau_Peierls_sc,ini%N_sl) allocate(prm%peierlsstress(prm%sum_N_sl,2)) - prm%peierlsstress(:,1) = prm%tau_Peierls_ed - prm%peierlsstress(:,2) = prm%tau_Peierls_sc + prm%peierlsstress(:,1) = math_expand(pl%get_as1dReal('tau_Peierls_ed', requiredSize=size(ini%N_sl)),ini%N_sl) + prm%peierlsstress(:,2) = math_expand(pl%get_as1dReal('tau_Peierls_sc', requiredSize=size(ini%N_sl)),ini%N_sl) - prm%rho_significant = pl%get_asReal('rho_significant') - prm%rho_min = pl%get_asReal('rho_min', 0.0_pREAL) - prm%C_CFL = pl%get_asReal('C_CFL',defaultVal=2.0_pREAL) + prm%rho_significant = pl%get_asReal('rho_significant') + prm%rho_min = pl%get_asReal('rho_min', 0.0_pREAL) + prm%C_CFL = pl%get_asReal('C_CFL',defaultVal=2.0_pREAL) - prm%V_at = pl%get_asReal('V_at') - prm%D_0 = pl%get_asReal('D_0') - prm%Q_cl = pl%get_asReal('Q_cl') - prm%f_F = pl%get_asReal('f_F') - prm%f_ed = pl%get_asReal('f_ed') - prm%w = pl%get_asReal('w') - prm%Q_sol = pl%get_asReal('Q_sol') - prm%f_sol = pl%get_asReal('f_sol') - prm%c_sol = pl%get_asReal('c_sol') + prm%V_at = pl%get_asReal('V_at') + prm%D_0 = pl%get_asReal('D_0') + prm%Q_cl = pl%get_asReal('Q_cl') + prm%f_F = pl%get_asReal('f_F') + prm%f_ed = pl%get_asReal('f_ed') + prm%w = pl%get_asReal('w') + prm%Q_sol = pl%get_asReal('Q_sol') + prm%f_sol = pl%get_asReal('f_sol') + prm%c_sol = pl%get_asReal('c_sol') - prm%p = pl%get_asReal('p_sl') - prm%q = pl%get_asReal('q_sl') - prm%B = pl%get_asReal('B') - prm%nu_a = pl%get_asReal('nu_a') + prm%p = pl%get_asReal('p_sl') + prm%q = pl%get_asReal('q_sl') + prm%B = pl%get_asReal('B') + prm%nu_a = pl%get_asReal('nu_a') ! ToDo: discuss logic - ini%sigma_rho_u = pl%get_asReal('sigma_rho_u') - ini%random_rho_u = pl%get_asReal('random_rho_u',defaultVal= 0.0_pREAL) + ini%sigma_rho_u = pl%get_asReal('sigma_rho_u') + ini%random_rho_u = pl%get_asReal('random_rho_u',defaultVal= 0.0_pREAL) if (pl%contains('random_rho_u')) & - ini%random_rho_u_binning = pl%get_asReal('random_rho_u_binning',defaultVal=0.0_pREAL) !ToDo: useful default? + ini%random_rho_u_binning = pl%get_asReal('random_rho_u_binning',defaultVal=0.0_pREAL) !ToDo: useful default? ! if (rhoSglRandom(instance) < 0.0_pREAL) & ! if (rhoSglRandomBinning(instance) <= 0.0_pREAL) & diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index b83c7a2d4..1fcc5eb69 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -139,7 +139,18 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) prm%sum_N_sl = sum(abs(N_sl)) slipActive: if (prm%sum_N_sl > 0) then - prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) + prm%dot_gamma_0_sl = pl%get_asReal('dot_gamma_0_sl') + prm%n_sl = pl%get_asReal('n_sl') + prm%a_sl = pl%get_asReal('a_sl') + prm%h_0_sl_sl = pl%get_asReal('h_0_sl-sl') + + xi_0_sl = math_expand(pl%get_as1dReal('xi_0_sl', requiredSize=size(N_sl)),N_sl) + prm%xi_inf_sl = math_expand(pl%get_as1dReal('xi_inf_sl', requiredSize=size(N_sl)),N_sl) + prm%h_int = math_expand(pl%get_as1dReal('h_int', requiredSize=size(N_sl), & + defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl) + + prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph)) + prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) if (phase_lattice(ph) == 'cI') then @@ -151,33 +162,21 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) prm%P_nS_pos = prm%P_sl prm%P_nS_neg = prm%P_sl end if - prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph)) - xi_0_sl = pl%get_as1dReal('xi_0_sl', requiredSize=size(N_sl)) - prm%xi_inf_sl = pl%get_as1dReal('xi_inf_sl', requiredSize=size(N_sl)) - prm%h_int = pl%get_as1dReal('h_int', requiredSize=size(N_sl), & - defaultVal=[(0.0_pREAL,i=1,size(N_sl))]) - - prm%dot_gamma_0_sl = pl%get_asReal('dot_gamma_0_sl') - prm%n_sl = pl%get_asReal('n_sl') - prm%a_sl = pl%get_asReal('a_sl') - prm%h_0_sl_sl = pl%get_asReal('h_0_sl-sl') - - ! expand: family => system - xi_0_sl = math_expand(xi_0_sl, N_sl) - prm%xi_inf_sl = math_expand(prm%xi_inf_sl,N_sl) - prm%h_int = math_expand(prm%h_int, N_sl) + prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) ! sanity checks - if ( prm%dot_gamma_0_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_sl' - if ( prm%a_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' a_sl' - if ( prm%n_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' n_sl' - if (any(xi_0_sl <= 0.0_pREAL)) extmsg = trim(extmsg)//' xi_0_sl' - if (any(prm%xi_inf_sl <= 0.0_pREAL)) extmsg = trim(extmsg)//' xi_inf_sl' + if ( prm%dot_gamma_0_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_sl' + if ( prm%a_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' a_sl' + if ( prm%n_sl <= 0.0_pREAL) extmsg = trim(extmsg)//' n_sl' + if (any(xi_0_sl <= 0.0_pREAL)) extmsg = trim(extmsg)//' xi_0_sl' + if (any(prm%xi_inf_sl <= 0.0_pREAL)) extmsg = trim(extmsg)//' xi_inf_sl' else slipActive xi_0_sl = emptyRealArray - allocate(prm%xi_inf_sl,prm%h_int,source=emptyRealArray) + allocate(prm%xi_inf_sl, & + prm%h_int, & + source=emptyRealArray) allocate(prm%h_sl_sl(0,0)) end if slipActive @@ -186,24 +185,22 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray) prm%sum_N_tw = sum(abs(N_tw)) twinActive: if (prm%sum_N_tw > 0) then - prm%systems_tw = lattice_labels_twin(N_tw,phase_lattice(ph)) - prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) - prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph)) + prm%c_1 = pl%get_asReal('c_1',defaultVal=0.0_pREAL) + prm%c_2 = pl%get_asReal('c_2',defaultVal=1.0_pREAL) + prm%c_3 = pl%get_asReal('c_3',defaultVal=0.0_pREAL) + prm%c_4 = pl%get_asReal('c_4',defaultVal=0.0_pREAL) + prm%dot_gamma_0_tw = pl%get_asReal('dot_gamma_0_tw') + prm%n_tw = pl%get_asReal('n_tw') + prm%f_sat_sl_tw = pl%get_asReal('f_sat_sl-tw') + prm%h_0_tw_tw = pl%get_asReal('h_0_tw-tw') + + xi_0_tw = math_expand(pl%get_as1dReal('xi_0_tw',requiredSize=size(N_tw)),N_tw) + prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) + prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph)) - xi_0_tw = pl%get_as1dReal('xi_0_tw',requiredSize=size(N_tw)) - - prm%c_1 = pl%get_asReal('c_1',defaultVal=0.0_pREAL) - prm%c_2 = pl%get_asReal('c_2',defaultVal=1.0_pREAL) - prm%c_3 = pl%get_asReal('c_3',defaultVal=0.0_pREAL) - prm%c_4 = pl%get_asReal('c_4',defaultVal=0.0_pREAL) - prm%dot_gamma_0_tw = pl%get_asReal('dot_gamma_0_tw') - prm%n_tw = pl%get_asReal('n_tw') - prm%f_sat_sl_tw = pl%get_asReal('f_sat_sl-tw') - prm%h_0_tw_tw = pl%get_asReal('h_0_tw-tw') - - ! expand: family => system - xi_0_tw = math_expand(xi_0_tw,N_tw) + prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) + prm%systems_tw = lattice_labels_twin(N_tw,phase_lattice(ph)) ! sanity checks if (prm%dot_gamma_0_tw <= 0.0_pREAL) extmsg = trim(extmsg)//' dot_gamma_0_tw'