From 37fa6c2e14ade5e50f5597bd442fabfd2b6dc91f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Apr 2012 17:28:08 +0000 Subject: [PATCH] merged code for python and spectral solver interfacing (shared most of it anyway). put functionality for getting header information (dimension, homogenization, resolution) in functions in mesh.f90 --- code/DAMASK_python_interface.f90 | 295 -------------- code/DAMASK_spectral.f90 | 107 +----- code/DAMASK_spectral_AL.f90 | 300 +++++++-------- code/DAMASK_spectral_interface.f90 | 238 ++++++------ code/damask.core.pyf | 4 +- code/mesh.f90 | 401 ++++++++++++-------- code/numerics.f90 | 3 +- code/{DAMASK_quit.f90 => spectral_quit.f90} | 2 +- 8 files changed, 534 insertions(+), 816 deletions(-) delete mode 100644 code/DAMASK_python_interface.f90 rename code/{DAMASK_quit.f90 => spectral_quit.f90} (98%) diff --git a/code/DAMASK_python_interface.f90 b/code/DAMASK_python_interface.f90 deleted file mode 100644 index 62b5430d4..000000000 --- a/code/DAMASK_python_interface.f90 +++ /dev/null @@ -1,295 +0,0 @@ -! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH -! -! This file is part of DAMASK, -! the Düsseldorf Advanced Material Simulation Kit. -! -! 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 . -! -!-------------------------------------------------------------------------------------------------- -!* $Id$ -!-------------------------------------------------------------------------------------------------- -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Interfacing between the spectral solver and the material subroutines provided -!! by DAMASK -!-------------------------------------------------------------------------------------------------- -module DAMASK_interface - - implicit none - private - character(len=64), parameter, public :: FEsolver = 'Spectral' !< Keyword for spectral solver - character(len=5), parameter, public :: inputFileExtension = '.geom' !< File extension for geometry description - character(len=4), parameter, public :: logFileExtension = '.log' !< Dummy variable as the spectral solver has no log - character(len=1024), private :: geometryParameter, & !< Interpretated parameter given at command line - loadcaseParameter !< Interpretated parameter given at command line - - public :: getSolverWorkingDirectoryName, & !< Interpretated parameter given at command line - getSolverJobName, & - getLoadCase, & - getLoadCaseName, & - getModelName, & - DAMASK_interface_init - private :: rectifyPath, & - makeRelativePath, & - getPathSep - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief initializes the solver by interpreting the command line arguments. Also writes -!! information on computation on screen -!-------------------------------------------------------------------------------------------------- - -subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn) - use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) - - implicit none - character(len=1024) :: hostName, & !< name of computer - userName !< name of user calling the executable - character(len=1024), intent(in) :: loadcaseParameterIn - character(len=1024), intent(in) :: geometryParameterIn - integer, dimension(8) :: dateAndTime ! type default integer - - geometryParameter = geometryParameterIn - loadcaseParameter = loadcaseParameterIn - - call date_and_time(values = dateAndTime) - - call GET_ENVIRONMENT_VARIABLE('HOST',hostName) - call GET_ENVIRONMENT_VARIABLE('USER',userName) - - write(6,*) - write(6,*) '<<<+- DAMASK_spectral_interface init -+>>>' - write(6,*) '$Id$' -#include "compilation_info.f90" - write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',& - dateAndTime(2),'/',& - dateAndTime(1) - write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',& - dateAndTime(6),':',& - dateAndTime(7) - write(6,*) 'Host Name: ', trim(hostName) - write(6,*) 'User Name: ', trim(userName) - write(6,*) 'Path Separator: ', getPathSep() - write(6,*) 'Geometry Parameter: ', trim(geometryParameter) - write(6,*) 'Loadcase Parameter: ', trim(loadcaseParameter) - -end subroutine DAMASK_interface_init - -!-------------------------------------------------------------------------------------------------- -!> @brief extract working directory from loadcase file possibly based on current working dir -!-------------------------------------------------------------------------------------------------- - character(len=1024) function getSolverWorkingDirectoryName() - - implicit none - character(len=1024) :: cwd - character :: pathSep - - pathSep = getPathSep() - - if (geometryParameter(1:1) == pathSep) then ! absolute path given as command line argument - getSolverWorkingDirectoryName = geometryParameter(1:scan(geometryParameter,pathSep,back=.true.)) - else - call getcwd(cwd) - getSolverWorkingDirectoryName = trim(cwd)//pathSep//geometryParameter(1:scan(geometryParameter,pathSep,back=.true.)) - endif - - getSolverWorkingDirectoryName = rectifyPath(getSolverWorkingDirectoryName) - -end function getSolverWorkingDirectoryName - - -!-------------------------------------------------------------------------------------------------- -!> @brief basename of geometry file from command line arguments -!-------------------------------------------------------------------------------------------------- -character(len=1024) function getSolverJobName() - - implicit none - getSolverJobName = trim(getModelName())//'_'//trim(getLoadCase()) - -end function getSolverJobName - - -!-------------------------------------------------------------------------------------------------- -!> @brief basename of geometry file from command line arguments -!-------------------------------------------------------------------------------------------------- -character(len=1024) function getModelName() - - use prec, only: pInt - - implicit none - character(len=1024) :: cwd - integer :: posExt,posSep - character :: pathSep - - pathSep = getPathSep() - posExt = scan(geometryParameter,'.',back=.true.) - posSep = scan(geometryParameter,pathSep,back=.true.) - - if (posExt <= posSep) posExt = len_trim(geometryParameter)+1 ! no extension present - getModelName = geometryParameter(1:posExt-1_pInt) ! path to geometry file (excl. extension) - - if (scan(getModelName,pathSep) /= 1) then ! relative path given as command line argument - call getcwd(cwd) - getModelName = rectifyPath(trim(cwd)//'/'//getModelName) - else - getModelName = rectifyPath(getModelName) - endif - - getModelName = makeRelativePath(getSolverWorkingDirectoryName(),& - getModelName) - -end function getModelName - - -!-------------------------------------------------------------------------------------------------- -!> @brief name of load case file exluding extension -!-------------------------------------------------------------------------------------------------- -character(len=1024) function getLoadCase() - - implicit none - integer :: posExt,posSep - character :: pathSep - - pathSep = getPathSep() - posExt = scan(loadcaseParameter,'.',back=.true.) - posSep = scan(loadcaseParameter,pathSep,back=.true.) - - if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1 ! no extension present - getLoadCase = loadcaseParameter(posSep+1:posExt-1) ! name of load case file exluding extension - -end function getLoadCase - - -!-------------------------------------------------------------------------------------------------- -!> @brief relative path of loadcase from command line arguments -!-------------------------------------------------------------------------------------------------- -character(len=1024) function getLoadcaseName() - - implicit none - character(len=1024) :: cwd - integer :: posExt = 0, posSep - character :: pathSep - - pathSep = getPathSep() - getLoadcaseName = loadcaseParameter - posExt = scan(getLoadcaseName,'.',back=.true.) - posSep = scan(getLoadcaseName,pathSep,back=.true.) - - if (posExt <= posSep) getLoadcaseName = trim(getLoadcaseName)//('.load') ! no extension present - if (scan(getLoadcaseName,pathSep) /= 1) then ! relative path given as command line argument - call getcwd(cwd) - getLoadcaseName = rectifyPath(trim(cwd)//pathSep//getLoadcaseName) - else - getLoadcaseName = rectifyPath(getLoadcaseName) - endif - - getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),& - getLoadcaseName) -end function getLoadcaseName - - -!-------------------------------------------------------------------------------------------------- -!> @brief remove ../ and ./ from path -!-------------------------------------------------------------------------------------------------- -function rectifyPath(path) - - implicit none - character(len=*) :: path - character(len=len_trim(path)) :: rectifyPath - character :: pathSep - integer :: i,j,k,l !no pInt - - pathSep = getPathSep() - - !remove ./ from path - l = len_trim(path) - rectifyPath = path - do i = l,3,-1 - if ( rectifyPath(i-1:i) == '.'//pathSep .and. rectifyPath(i-2:i-2) /= '.' ) & - rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' - enddo - - !remove ../ and corresponding directory from rectifyPath - l = len_trim(rectifyPath) - i = index(rectifyPath(i:l),'..'//pathSep) - j = 0 - do while (i > j) - j = scan(rectifyPath(1:i-2),pathSep,back=.true.) - rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j) - if (rectifyPath(j+1:j+1) == pathSep) 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) = ' ' - endif - i = j+index(rectifyPath(j+1:l),'..'//pathSep) - enddo - if(len_trim(rectifyPath) == 0) rectifyPath = pathSep - -end function rectifyPath - - -!-------------------------------------------------------------------------------------------------- -!> @brief relative path from absolute a to absolute b -!-------------------------------------------------------------------------------------------------- -character(len=1024) function makeRelativePath(a,b) - - implicit none - character (len=*) :: a,b - character :: pathSep - integer :: i,posLastCommonSlash,remainingSlashes !no pInt - - pathSep = getPathSep() - posLastCommonSlash = 0 - remainingSlashes = 0 - - do i = 1, min(1024,len_trim(a),len_trim(b)) - if (a(i:i) /= b(i:i)) exit - if (a(i:i) == pathSep) posLastCommonSlash = i - enddo - do i = posLastCommonSlash+1,len_trim(a) - if (a(i:i) == pathSep) remainingSlashes = remainingSlashes + 1 - enddo - makeRelativePath = repeat('..'//pathSep,remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) - -end function makeRelativePath - - -!-------------------------------------------------------------------------------------------------- -!> @brief counting / and \ in $PATH System variable the character occuring more often is assumed -!! to be the path separator -!-------------------------------------------------------------------------------------------------- -character function getPathSep() - - use prec, only: pInt - - implicit none - character(len=2048) path - integer(pInt) :: backslash = 0_pInt, slash = 0_pInt - integer :: i - - call get_environment_variable('PATH',path) - do i=1, len(trim(path)) - if (path(i:i)=='/') slash = slash + 1_pInt - if (path(i:i)=='\') backslash = backslash + 1_pInt - enddo - - if (backslash>slash) then - getPathSep = '\' - else - getPathSep = '/' - endif - -end function - -end module diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 07cb36d16..86461f32b 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -34,8 +34,9 @@ ! ! MPI fuer Eisenforschung, Duesseldorf -program DAMASK_spectral +#include "spectral_quit.f90" +program DAMASK_spectral use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use DAMASK_interface, only: & @@ -75,6 +76,11 @@ program DAMASK_spectral use math + use mesh, only : & + mesh_spectral_getResolution, & + mesh_spectral_getDimension, & + mesh_spectral_getHomogenization + use CPFEM, only: & CPFEM_general, & CPFEM_initAll @@ -86,6 +92,7 @@ program DAMASK_spectral use numerics, only: & err_div_tol, & err_stress_tolrel, & + err_stress_tolabs, & rotation_tol, & itmax,& itmin, & @@ -102,11 +109,7 @@ program DAMASK_spectral !$ use OMP_LIB ! the openMP function library - !################################################################################################## -! variable declaration -!################################################################################################## implicit none - !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file real(pReal), dimension(9) :: & @@ -121,7 +124,6 @@ program DAMASK_spectral integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing integer(pInt) :: & - headerLength, & N_l = 0_pInt, & N_t = 0_pInt, & N_n = 0_pInt, & @@ -132,12 +134,6 @@ program DAMASK_spectral character(len=1024) :: & line, & - keyword - - logical :: & - gotResolution = .false.,& - gotDimension = .false.,& - gotHomogenization = .false. type bc_type real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient @@ -252,7 +248,6 @@ program DAMASK_spectral !################################################################################################## ! reading of information from load case file and geometry file !################################################################################################## - !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS open (6, encoding='UTF-8') call DAMASK_interface_init write(6,'(a)') '' @@ -360,72 +355,16 @@ program DAMASK_spectral end select enddo; enddo 101 close(myUnit) -if (sum(bc(1:N_Loadcases)%incs)>9000_pInt) stop 'to many incs' !discuss with Philip, stop code trouble. suggesting warning !-------------------------------------------------------------------------------------------------- ToDo: if temperature at CPFEM is treated properly, move this up immediately after interface init ! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry) call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt) - if (update_gamma .and. .not. memory_efficient) call IO_error(error_ID = 847_pInt) - + !-------------------------------------------------------------------------------------------------- -! read header of geom file to get size information. complete geom file is intepretated by mesh.f90 - call IO_open_file(myUnit,trim(getModelName())//InputFileExtension) - rewind(myUnit) - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,2_pInt) - keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) - if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt - else - call IO_error(error_ID=842_pInt) - endif - rewind(myUnit) - do i = 1_pInt, headerLength - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,maxNchunksGeom) - select case ( IO_lc(IO_StringValue(line,positions,1)) ) - case ('dimension') - gotDimension = .true. - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,positions,j))) - case('x') - geomdim(1) = IO_floatValue(line,positions,j+1_pInt) - case('y') - geomdim(2) = IO_floatValue(line,positions,j+1_pInt) - case('z') - geomdim(3) = IO_floatValue(line,positions,j+1_pInt) - end select - enddo - case ('homogenization') - gotHomogenization = .true. - homog = IO_intValue(line,positions,2_pInt) - case ('resolution') - gotResolution = .true. - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,positions,j))) - case('a') - res(1) = IO_intValue(line,positions,j+1_pInt) - case('b') - res(2) = IO_intValue(line,positions,j+1_pInt) - case('c') - res(3) = IO_intValue(line,positions,j+1_pInt) - end select - enddo - end select - enddo - close(myUnit) - -!-------------------------------------------------------------------------------------------------- -! sanity checks of geometry parameters - if (.not.(gotDimension .and. gotHomogenization .and. gotResolution))& - call IO_error(error_ID = 845_pInt) - if (any(geomdim<=0.0_pReal)) call IO_error(error_ID = 802_pInt) - if(mod(res(1),2_pInt)/=0_pInt .or.& - mod(res(2),2_pInt)/=0_pInt .or.& - (mod(res(3),2_pInt)/=0_pInt .and. res(3)/= 1_pInt)) call IO_error(error_ID = 803_pInt) - -!-------------------------------------------------------------------------------------------------- -! variables derived from resolution +! get resolution, dimension, homogenization and variables derived from resolution + res = mesh_spectral_getResolution(myUnit) + geomdim = mesh_spectral_getDimension(myUnit) + homog = mesh_spectral_getHomogenization(myUnit) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) Npoints = res(1)*res(2)*res(3) wgt = 1.0_pReal/real(Npoints, pReal) @@ -519,13 +458,11 @@ if (sum(bc(1:N_Loadcases)%incs)>9000_pInt) stop 'to many incs' !discuss with Phi !-------------------------------------------------------------------------------------------------- ! general initialization of fftw (see manual on fftw.org for more details) if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C -#ifdef _OPENMP - if(DAMASK_NumThreadsInt > 0_pInt) then - ierr = fftw_init_threads() - if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt) - call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) - endif -#endif +!$ if(DAMASK_NumThreadsInt > 0_pInt) then +!$ ierr = fftw_init_threads() +!$ if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt) +!$ call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) +!$ endif call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation !-------------------------------------------------------------------------------------------------- @@ -880,7 +817,7 @@ C_ref = C * wgt ! stress BC handling if(size_reduced > 0_pInt) then ! calculate stress BC if applied err_stress = maxval(abs(mask_stress * (P_av - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable) - err_stress_tol = maxval(abs(P_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent + err_stress_tol = min(maxval(abs(P_av)) * err_stress_tolrel,err_stress_tolabs) ! don't use any tensor norm for the relative criterion because the comparison should be coherent write(6,'(a)') '' write(6,'(a)') '... correcting deformation gradient to fulfill BCs ...............' write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, & @@ -1124,7 +1061,5 @@ C_ref = C * wgt call fftw_destroy_plan(plan_scalarField_forth) call fftw_destroy_plan(plan_scalarField_back) endif - call quit(0_pInt) -end program DAMASK_spectral - -#include "DAMASK_quit.f90" \ No newline at end of file + call quit(1_pInt) +end program DAMASK_spectral \ No newline at end of file diff --git a/code/DAMASK_spectral_AL.f90 b/code/DAMASK_spectral_AL.f90 index 0a42175fe..4adbe9ab9 100644 --- a/code/DAMASK_spectral_AL.f90 +++ b/code/DAMASK_spectral_AL.f90 @@ -33,32 +33,82 @@ ! R. Lebensohn ! ! MPI fuer Eisenforschung, Duesseldorf -!################################################################################################## -! used modules -!################################################################################################## -program DAMASK_spectral_AL +#include "spectral_quit.f90" +program DAMASK_spectral_AL use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) - use DAMASK_interface - use prec, only: pInt, pReal, DAMASK_NaN - use IO - use debug, only: debug_spectral, & - debug_levelBasic, & - debug_spectralRestart, & - debug_spectralFFTW - use math - use CPFEM, only: CPFEM_general, CPFEM_initAll - use FEsolving, only: restartWrite, restartInc - use numerics, only: err_div_tol, err_stress_tolrel, rotation_tol, itmax, itmin, & - memory_efficient, update_gamma, DAMASK_NumThreadsInt, & - fftw_planner_flag, fftw_timelimit - use homogenization, only: materialpoint_sizeResults, materialpoint_results - !$ use OMP_LIB ! the openMP function library -!################################################################################################## -! variable declaration -!################################################################################################## - implicit none + use DAMASK_interface, only: & + DAMASK_interface_init, & + getLoadcaseName, & + getSolverWorkingDirectoryName, & + getSolverJobName, & + getModelName, & + inputFileExtension + + use prec, only: & + pInt, & + pReal, & + DAMASK_NaN + + use IO, only: & + IO_isBlank, & + IO_open_file, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_error, & + IO_lc, & + IO_read_jobBinaryFile, & + IO_write_jobBinaryFile + + use debug, only: & + debug_what, & + debug_spectral, & + debug_levelBasic, & + debug_spectralDivergence, & + debug_spectralRestart, & + debug_spectralFFTW, & + debug_reset, & + debug_info + + use math + + use mesh, only : & + mesh_spectral_getResolution, & + mesh_spectral_getDimension, & + mesh_spectral_getHomogenization + + use CPFEM, only: & + CPFEM_general, & + CPFEM_initAll + + use FEsolving, only: & + restartWrite, & + restartInc + + use numerics, only: & + err_div_tol, & + err_stress_tolrel, & + err_stress_tolabs, & + rotation_tol, & + itmax,& + itmin, & + memory_efficient, & + update_gamma, & + divergence_correction, & + DAMASK_NumThreadsInt, & + fftw_planner_flag, & + fftw_timelimit + + use homogenization, only: & + materialpoint_sizeResults, & + materialpoint_results + + !$ use OMP_LIB ! the openMP function library + + implicit none !-------------------------------------------------------------------------------------------------- ! variables to read from load case and geom file real(pReal), dimension(9) :: temp_valueVector ! stores information temporarily from loadcase file @@ -69,15 +119,12 @@ program DAMASK_spectral_AL maxNchunksGeom = 7_pInt, & ! 4 identifiers, 3 values myUnit = 234_pInt integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing - integer(pInt) :: headerLength,& - N_l = 0_pInt,& - N_t = 0_pInt,& - N_n = 0_pInt,& - N_Fdot = 0_pInt - character(len=1024) :: path, line, keyword - logical :: gotResolution = .false.,& - gotDimension = .false.,& - gotHomogenization = .false. + integer(pInt) :: & + N_l = 0_pInt,& + N_t = 0_pInt,& + N_n = 0_pInt,& + N_Fdot = 0_pInt + character(len=1024) :: line !-------------------------------------------------------------------------------------------------- ! variable storing information from load case file @@ -161,21 +208,19 @@ program DAMASK_spectral_AL !################################################################################################## ! reading of information from load case file and geometry file !################################################################################################## - !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS open (6, encoding='UTF-8') call DAMASK_interface_init - print '(a)', '' - print '(a)', ' <<<+- DAMASK_spectral_AL init -+>>>' - print '(a)', ' $Id$' + write(6,'(a)') '' + write(6,'(a)') ' <<<+- DAMASK_spectral_AL init -+>>>' + write(6,'(a)') ' $Id$' #include "compilation_info.f90" - print '(a,a)', ' Working Directory: ',trim(getSolverWorkingDirectoryName()) - print '(a,a)', ' Solver Job Name: ',trim(getSolverJobName()) - print '(a)', '' + write(6,'(a)') ' Working Directory: ',trim(getSolverWorkingDirectoryName()) + write(6,'(a)') ' Solver Job Name: ',trim(getSolverJobName()) + write(6,'(a)') '' !-------------------------------------------------------------------------------------------------- ! reading the load case file and allocate data structure containing load cases - path = getLoadcaseName() - call IO_open_file(myUnit,path) + call IO_open_file(myUnit,getLoadcaseName()) rewind(myUnit) do read(myUnit,'(a1024)',END = 100) line @@ -197,7 +242,7 @@ program DAMASK_spectral_AL 100 N_Loadcases = N_n if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check - call IO_error(error_ID=837_pInt,ext_msg = trim(path)) ! error message for incomplete loadcase + call IO_error(error_ID=837_pInt,ext_msg = trim(getLoadcaseName())) ! error message for incomplete loadcase allocate (bc(N_Loadcases)) !-------------------------------------------------------------------------------------------------- @@ -274,88 +319,31 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ToDo: if temperature at CPFEM is treated properly, move this up immediately after interface init ! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry) call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt) - if (update_gamma .and. .not. memory_efficient) call IO_error(error_ID = 847_pInt) - -!-------------------------------------------------------------------------------------------------- -! read header of geom file to get size information. complete geom file is intepretated by mesh.f90 - path = getModelName() - call IO_open_file(myUnit,trim(path)//InputFileExtension) - rewind(myUnit) - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,2_pInt) - keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) - if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt - else - call IO_error(error_ID=842_pInt) - endif - rewind(myUnit) - do i = 1_pInt, headerLength - read(myUnit,'(a1024)') line - positions = IO_stringPos(line,maxNchunksGeom) - select case ( IO_lc(IO_StringValue(line,positions,1)) ) - case ('dimension') - gotDimension = .true. - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,positions,j))) - case('x') - geomdim(1) = IO_floatValue(line,positions,j+1_pInt) - case('y') - geomdim(2) = IO_floatValue(line,positions,j+1_pInt) - case('z') - geomdim(3) = IO_floatValue(line,positions,j+1_pInt) - end select - enddo - case ('homogenization') - gotHomogenization = .true. - homog = IO_intValue(line,positions,2_pInt) - case ('resolution') - gotResolution = .true. - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,positions,j))) - case('a') - res(1) = IO_intValue(line,positions,j+1_pInt) - case('b') - res(2) = IO_intValue(line,positions,j+1_pInt) - case('c') - res(3) = IO_intValue(line,positions,j+1_pInt) - end select - enddo - end select - enddo - close(myUnit) - !-------------------------------------------------------------------------------------------------- -! sanity checks of geometry parameters - if (.not.(gotDimension .and. gotHomogenization .and. gotResolution))& - call IO_error(error_ID = 845_pInt) - if (any(geomdim<=0.0_pReal)) call IO_error(error_ID = 802_pInt) - if(mod(res(1),2_pInt)/=0_pInt .or.& - mod(res(2),2_pInt)/=0_pInt .or.& - (mod(res(3),2_pInt)/=0_pInt .and. res(3)/= 1_pInt)) call IO_error(error_ID = 803_pInt) - -!-------------------------------------------------------------------------------------------------- -! variables derived from resolution +! get resolution, dimension, homogenization and variables derived from resolution + res = mesh_spectral_getResolution(myUnit) + geomdim = mesh_spectral_getDimension(myUnit) + homog = mesh_spectral_getHomogenization(myUnit) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) Npoints = res(1)*res(2)*res(3) wgt = 1.0_pReal/real(Npoints, pReal) - + !-------------------------------------------------------------------------------------------------- ! output of geometry - print '(a)', '' - print '(a)', '#############################################################' - print '(a)', 'DAMASK spectral_AL:' - print '(a)', 'The AL spectral method boundary value problem solver for' - print '(a)', 'the Duesseldorf Advanced Material Simulation Kit' - print '(a)', '#############################################################' - print '(a,a)', 'geometry file: ',trim(path)//'.geom' - print '(a)', '=============================================================' - print '(a,3(i12 ))','resolution a b c:', res - print '(a,3(f12.5))','dimension x y z:', geomdim - print '(a,i5)','homogenization: ',homog - print '(a)', '#############################################################' - print '(a,a)', 'loadcase file: ',trim(getLoadcaseName()) + write(6,'(a)') '' + write(6,'(a)') '#############################################################' + write(6,'(a)') 'DAMASK spectral_AL:' + write(6,'(a)') 'The AL spectral method boundary value problem solver for' + write(6,'(a)') 'the Duesseldorf Advanced Material Simulation Kit' + write(6,'(a)') '#############################################################' + write(6,'(a)') 'geometry file: ',trim(getModelName())//InputFileExtension + write(6,'(a)') '=============================================================' + write(6,'(a,3(i12 ))') 'resolution a b c:', res + write(6,'(a,3(f12.5))') 'dimension x y z:', geomdim + write(6,'(a,i5)') 'homogenization: ',homog + write(6,'(a)') '#############################################################' + write(6,'(a)') 'loadcase file: ',trim(getLoadcaseName()) !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case @@ -364,32 +352,32 @@ program DAMASK_spectral_AL do loadcase = 1_pInt, N_Loadcases write (loadcase_string, '(i6)' ) loadcase - print '(a)', '=============================================================' - print '(a,i6)', 'loadcase: ', loadcase + write(6,'(a)') '=============================================================' + write(6,'(a,i6)') 'loadcase: ', loadcase - if (.not. bc(loadcase)%followFormerTrajectory) print '(a)', 'drop guessing along trajectory' + if (.not. bc(loadcase)%followFormerTrajectory) write(6,'(a)') 'drop guessing along trajectory' if (bc(loadcase)%velGradApplied) then do j = 1_pInt, 3_pInt if (any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .true.) .and. & any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined enddo - print '(a)','velocity gradient:' + write(6,'(a)') 'velocity gradient:' else - print '(a)','deformation gradient rate:' + write(6,'(a)') 'deformation gradient rate:' endif - write (*,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),& + write (6,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),& reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskDeformation)) - write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'stress / GPa:',& + write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'stress / GPa:',& 1e-9_pReal*merge(math_transpose33(bc(loadcase)%P),& reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskStress)) if (any(bc(loadcase)%rotation /= math_I3)) & write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',& math_transpose33(bc(loadcase)%rotation) - print '(a,f12.6)','temperature:',bc(loadcase)%temperature - print '(a,f12.6)','time: ',bc(loadcase)%time - print '(a,i5)' ,'increments: ',bc(loadcase)%incs - print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency - print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency + write(6,'(a,f12.6)') 'temperature:',bc(loadcase)%temperature + write(6,'(a,f12.6)') 'time: ',bc(loadcase)%time + write(6,'(a,i5)') 'increments: ',bc(loadcase)%incs + write(6,'(a,i5)') 'output frequency: ',bc(loadcase)%outputfrequency + write(6,'(a,i5)') 'restart frequency: ',bc(loadcase)%restartfrequency if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 831_pInt ! exclusive or masking only if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & @@ -434,13 +422,11 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! general initialization of fftw (see manual on fftw.org for more details) if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C -#ifdef _OPENMP - if(DAMASK_NumThreadsInt > 0_pInt) then - ierr = fftw_init_threads() - if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt) - call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) - endif -#endif +!$ if(DAMASK_NumThreadsInt > 0_pInt) then +!$ ierr = fftw_init_threads() +!$ if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt) +!$ call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) +!$ endif call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation !-------------------------------------------------------------------------------------------------- @@ -456,7 +442,7 @@ program DAMASK_spectral_AL 1, res(3)*res(2)* res1_red,& F_real,[ res(3),res(2) ,res(1)+2_pInt],& 1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) - if (debugGeneral) print '(a)' , 'FFTW initialized' + if (debugGeneral) write(6,'(a)') 'FFTW initialized' !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) and remove the given highest frequencies @@ -514,7 +500,7 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! possible restore deformation gradient from saved state if (restartInc > 1_pInt) then ! using old values from file - if (debugRestart) print '(a,i6,a)' , 'Reading values of increment ',& + if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& restartInc - 1_pInt,' from file' call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& trim(getSolverJobName()),size(F_star)) @@ -548,7 +534,7 @@ program DAMASK_spectral_AL write(538) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc write(538) 'eoh' ! end of header write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results - if (debugGeneral) print '(a)' , 'Header of result file written out' + if (debugGeneral) write(6,'(a)') 'Header of result file written out' !################################################################################################## ! Loop over loadcases defined in the loadcase file @@ -677,8 +663,8 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! report begin of new increment - print '(a)', '##################################################################' - print '(A,I5.5,A,es12.5)', 'Increment ', totalIncsCounter, ' Time ',time + write(6,'(a)') '##################################################################' + write(6,'(A,I5.5,A,es12.5)') 'Increment ', totalIncsCounter, ' Time ',time guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase CPFEM_mode = 1_pInt ! winding forward @@ -696,9 +682,9 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! report begin of new iteration - print '(a)', '' - print '(a)', '==================================================================' - print '(5(a,i6.6))', 'Loadcase ',loadcase,' Increment ',inc,'/',bc(loadcase)%incs,& + write(6,'(a)') '' + write(6,'(a)') '==================================================================' + write(6,'(5(a,i6.6))') 'Loadcase ',loadcase,' Increment ',inc,'/',bc(loadcase)%incs,& ' @ Iteration ',iter,'/',itmax !-------------------------------------------------------------------------------------------------- @@ -716,7 +702,7 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! doing Fourier transform - print '(a)', '... spectral method ...............................................' + write(6,'(a)') '... spectral method ...............................................' lambda_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = lambda(1:res(1),1:res(2),1:res(3),1:3,1:3) call fftw_execute_dft_r2c(plan_lambda,lambda_real,lambda_fourier) lambda_fourier( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& @@ -751,10 +737,6 @@ program DAMASK_spectral_AL enddo; enddo err_div_RMS = sqrt(err_div_RMS)*wgt - ! if (err_div < err_div_RMS/pstress_av_L2 .and. guessmax<0) then - ! print*, 'increasing div, stopping calc' - ! iter = huge(1_pInt) - ! endif err_div = err_div_RMS/pstress_av_L2 !-------------------------------------------------------------------------------------------------- ! using gamma operator to update F @@ -793,7 +775,7 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! if(callCPFEM) then - print '(a)', '... calling CPFEM to update P(F*) and F*.........................' + write(6,'(a)') '... calling CPFEM to update P(F*) and F*.........................' F_star_lastIter = F_star ielem = 0_pInt do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) @@ -819,7 +801,7 @@ program DAMASK_spectral_AL enddo; enddo; enddo else guesses = guesses +1_pInt - print*, '... linear approximation for P(F*) and F* ', guesses, ' of ', guessmax + write(6,'(a)')' ... linear approximation for P(F*) and F* ', guesses, ' of ', guessmax do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) temp33_Real = lambda(i,j,k,1:3,1:3) - (P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3),& F_star(i,j,k,1:3,1:3) -F_star_lastIter(i,j,k,1:3,1:3)))& @@ -830,7 +812,7 @@ program DAMASK_spectral_AL enddo; enddo; enddo endif - print '(a)', '... update λ..........................' + write(6,'(a)') '... update λ..........................' err_f = 0.0_pReal err_f_point = 0.0_pReal @@ -879,7 +861,7 @@ program DAMASK_spectral_AL write(6,'(a,es14.7)') 'max abs err F', err_f_point write(6,'(a,es14.7)') 'max abs err P', err_p_point err_crit = max(err_p/1e-3, err_f/1e-4,err_div/err_div_tol,err_stress/err_stress_tol) - print*, 'critical error', err_crit + write(6,'(a)') 'critical error', err_crit if (.not. callCPFEM) then if(err_crit < 1.0_pReal .or. guesses >= guessmax) callCPFEM = .true. @@ -920,13 +902,11 @@ program DAMASK_spectral_AL deallocate(c_reduced) deallocate(s_reduced) enddo ! end looping over loadcases - print '(a)', '' - print '(a)', '##################################################################' - print '(i6.6,a,i6.6,a)', notConvergedCounter, ' out of ', & - notConvergedCounter + convergedCounter, ' increments did not converge!' + write(6,'(a)') '' + write(6,'(a)') '##################################################################' + write(6,'(i6.6,a,i6.6,a)') notConvergedCounter, ' out of ', & + notConvergedCounter + convergedCounter, ' increments did not converge!' close(538) call fftw_destroy_plan(plan_lambda); call fftw_destroy_plan(plan_correction) call quit(1_pInt) -end program DAMASK_spectral_AL - -#include "DAMASK_quit.f90" \ No newline at end of file +end program DAMASK_spectral_AL \ No newline at end of file diff --git a/code/DAMASK_spectral_interface.f90 b/code/DAMASK_spectral_interface.f90 index 9979c7bed..292e43164 100644 --- a/code/DAMASK_spectral_interface.f90 +++ b/code/DAMASK_spectral_interface.f90 @@ -38,7 +38,7 @@ module DAMASK_interface getLoadCase, & getLoadCaseName, & getModelName, & - DAMASK_interface_init + DAMASK_interface_init private :: rectifyPath, & makeRelativePath, & getPathSep @@ -49,124 +49,140 @@ contains !> @brief initializes the solver by interpreting the command line arguments. Also writes !! information on computation on screen !-------------------------------------------------------------------------------------------------- -subroutine DAMASK_interface_init +subroutine DAMASK_interface_init(loadcaseParameterIn,geometryParameterIn) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pInt implicit none - character(len=1024) :: commandLine, & !< command line call as string - hostName, & !< name of computer - userName !< name of user calling the executable - integer :: i, & - start ,& - length - integer, dimension(8) :: dateAndTime ! type default integer - - call get_command(commandLine) - call date_and_time(values = dateAndTime) - do i = 1,len(commandLine) ! remove capitals - if(64 0 .or. index(commandLine,' --help ',.true.) > 0) then ! search for ' -h ' or '--help' - write(6,'(a)') '$Id$' -#include "compilation_info.f90" - write(6,'(a)') '#############################################################' - write(6,'(a)') 'DAMASK spectral:' - write(6,'(a)') 'The spectral method boundary value problem solver for' - write(6,'(a)') 'the Duesseldorf Advanced Material Simulation Kit' - write(6,'(a)') '#############################################################' - write(6,'(a)') 'Valid command line switches:' - write(6,'(a)') ' --geom (-g, --geometry)' - write(6,'(a)') ' --load (-l, --loadcase)' - write(6,'(a)') ' --restart (-r)' - write(6,'(a)') ' --help (-h)' - write(6,'(a)') ' ' - write(6,'(a)') 'Mandatory Arguments:' - write(6,'(a)') ' --load PathToLoadFile/NameOfLoadFile.load' - write(6,'(a)') ' "PathToGeomFile" will be the working directory.' - write(6,'(a)') ' Make sure the file "material.config" exists in the working' - write(6,'(a)') ' directory' - write(6,'(a)') ' For further configuration place "numerics.config"' - write(6,'(a)') ' and "numerics.config" in that directory.' - write(6,'(a)') ' ' - write(6,'(a)') ' --geom PathToGeomFile/NameOfGeom.geom' - write(6,'(a)') ' ' - write(6,'(a)') 'Optional Argument:' - write(6,'(a)') ' --restart XX' - write(6,'(a)') ' Reads in total increment No. XX-1 and continous to' - write(6,'(a)') ' calculate total increment No. XX.' - write(6,'(a)') ' Attention: Overwrites existing results file ' - write(6,'(a)') ' "NameOfGeom_NameOfLoadFile_spectralOut".' - write(6,'(a)') ' Works only if the restart information for total increment' - write(6,'(a)') ' No. XX-1 is available in the working directory.' - write(6,'(a)') 'Help:' - write(6,'(a)') ' --help' - write(6,'(a)') ' Prints this message and exits' - write(6,'(a)') ' ' - call quit(0_pInt) - endif - if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) & ! check for correct number of given arguments (no --help) - stop 'Wrong Nr. of Arguments. Run DAMASK_spectral.exe --help' ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available - start = index(commandLine,'-g',.true.) + 3 ! search for '-g' and jump to first char of geometry - if (index(commandLine,'--geom',.true.)>0) then ! if '--geom' is found, use that (contains '-g') - start = index(commandLine,'--geom',.true.) + 7 - endif - if (index(commandLine,'--geometry',.true.)>0) then ! again, now searching for --geometry' - start = index(commandLine,'--geometry',.true.) + 11 - endif - if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available - write(6,'(a)') 'No Geometry specified' - call quit(9999) - endif - length = index(commandLine(start:len(commandLine)),' ',.false.) - - call get_command(commandLine) ! may contain capitals - geometryParameter = '' ! should be empty - geometryParameter(1:length)=commandLine(start:start+length) - - do i=1,len(commandLine) ! remove capitals - if(640) then ! if '--load' is found, use that (contains '-l') - start = index(commandLine,'--load',.true.) + 7 - endif - if (index(commandLine,'--loadcase',.true.)>0) then ! again, now searching for --loadcase' - start = index(commandLine,'--loadcase',.true.) + 11 - endif - if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available - write(6,'(a)') 'No Loadcase specified' - call quit(9999) - endif - length = index(commandLine(start:len(commandLine)),' ',.false.) - - call get_command(commandLine) ! may contain capitals - loadcaseParameter = '' ! should be empty - loadcaseParameter(1:length)=commandLine(start:start+length) - - do i=1,len(commandLine) ! remove capitals - if(640) then ! if '--restart' is found, use that (contains '-l') - start = index(commandLine,'--restart',.true.) + 7 - endif - length = index(commandLine(start:len(commandLine)),' ',.false.) - - call get_command(commandLine) ! may contain capitals - call GET_ENVIRONMENT_VARIABLE('HOST',hostName) - call GET_ENVIRONMENT_VARIABLE('USER',userName) - + character(len=1024), optional, intent(in) :: & + loadcaseParameterIn, & + geometryParameterIn + character(len=1024) :: & + commandLine, & !< command line call as string + hostName, & !< name of computer + userName !< name of user calling the executable + integer :: & + i, & + start ,& + length + integer, dimension(8) :: & + dateAndTime ! type default integer write(6,*) write(6,*) '<<<+- DAMASK_spectral_interface init -+>>>' write(6,*) '$Id$' #include "compilation_info.f90" + if ( present(loadcaseParameterIn) .and. present(geometryParameterIn)) then ! both mandatory parameters given in function call + geometryParameter = geometryParameterIn + loadcaseParameter = loadcaseParameterIn + commandLine='n/a' + start = 3_pInt + else if ( .not.( present(loadcaseParameterIn) .and. present(geometryParameterIn))) then ! none parameters given in function call, trying to get them from comman line + call get_command(commandLine) + call date_and_time(values = dateAndTime) + do i = 1,len(commandLine) ! remove capitals + if(64 0 .or. index(commandLine,' --help ',.true.) > 0) then ! search for ' -h ' or '--help' + write(6,'(a)') '#############################################################' + write(6,'(a)') 'DAMASK spectral:' + write(6,'(a)') 'The spectral method boundary value problem solver for' + write(6,'(a)') 'the Duesseldorf Advanced Material Simulation Kit' + write(6,'(a)') '#############################################################' + write(6,'(a)') 'Valid command line switches:' + write(6,'(a)') ' --geom (-g, --geometry)' + write(6,'(a)') ' --load (-l, --loadcase)' + write(6,'(a)') ' --restart (-r)' + write(6,'(a)') ' --help (-h)' + write(6,'(a)') ' ' + write(6,'(a)') 'Mandatory Arguments:' + write(6,'(a)') ' --load PathToLoadFile/NameOfLoadFile.load' + write(6,'(a)') ' "PathToGeomFile" will be the working directory.' + write(6,'(a)') ' Make sure the file "material.config" exists in the working' + write(6,'(a)') ' directory' + write(6,'(a)') ' For further configuration place "numerics.config"' + write(6,'(a)') ' and "numerics.config" in that directory.' + write(6,'(a)') ' ' + write(6,'(a)') ' --geom PathToGeomFile/NameOfGeom.geom' + write(6,'(a)') ' ' + write(6,'(a)') 'Optional Argument:' + write(6,'(a)') ' --restart XX' + write(6,'(a)') ' Reads in total increment No. XX-1 and continous to' + write(6,'(a)') ' calculate total increment No. XX.' + write(6,'(a)') ' Attention: Overwrites existing results file ' + write(6,'(a)') ' "NameOfGeom_NameOfLoadFile_spectralOut".' + write(6,'(a)') ' Works only if the restart information for total increment' + write(6,'(a)') ' No. XX-1 is available in the working directory.' + write(6,'(a)') 'Help:' + write(6,'(a)') ' --help' + write(6,'(a)') ' Prints this message and exits' + write(6,'(a)') ' ' + call quit(1_pInt) ! normal Termination + endif + if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) then ! check for correct number of given arguments (no --help) + write(6,'(a)') 'Wrong Nr. of Arguments. Run DAMASK_spectral.exe --help' ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available + call quit(100_pInt) ! abnormal termination + endif + start = index(commandLine,'-g',.true.) + 3 ! search for '-g' and jump to first char of geometry + if (index(commandLine,'--geom',.true.)>0) then ! if '--geom' is found, use that (contains '-g') + start = index(commandLine,'--geom',.true.) + 7 + endif + if (index(commandLine,'--geometry',.true.)>0) then ! again, now searching for --geometry' + start = index(commandLine,'--geometry',.true.) + 11 + endif + if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available + write(6,'(a)') 'No Geometry specified' + call quit(100_pInt) ! abnormal termination + endif + length = index(commandLine(start:len(commandLine)),' ',.false.) + + call get_command(commandLine) ! may contain capitals + geometryParameter = '' ! should be empty + geometryParameter(1:length)=commandLine(start:start+length) + + do i=1,len(commandLine) ! remove capitals + if(640) then ! if '--load' is found, use that (contains '-l') + start = index(commandLine,'--load',.true.) + 7 + endif + if (index(commandLine,'--loadcase',.true.)>0) then ! again, now searching for --loadcase' + start = index(commandLine,'--loadcase',.true.) + 11 + endif + if(start==3_pInt) then ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available + write(6,'(a)') 'No Loadcase specified' + call quit(100_pInt) ! abnormal termination + endif + length = index(commandLine(start:len(commandLine)),' ',.false.) + + call get_command(commandLine) ! may contain capitals + loadcaseParameter = '' ! should be empty + loadcaseParameter(1:length)=commandLine(start:start+length) + + do i=1,len(commandLine) ! remove capitals + if(640) then ! if '--restart' is found, use that (contains '-l') + start = index(commandLine,'--restart',.true.) + 7 + endif + length = index(commandLine(start:len(commandLine)),' ',.false.) + + call get_command(commandLine) ! may contain capitals + else + write(6,'(a)') 'Wrong Nr. of Arguments!' ! Function call with wrong No of arguments + call quit(100_pInt) + endif + + call GET_ENVIRONMENT_VARIABLE('HOST',hostName) + call GET_ENVIRONMENT_VARIABLE('USER',userName) + write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',& dateAndTime(2),'/',& dateAndTime(1) diff --git a/code/damask.core.pyf b/code/damask.core.pyf index d2c48d0ba..b6be470ff 100644 --- a/code/damask.core.pyf +++ b/code/damask.core.pyf @@ -15,9 +15,9 @@ python module core ! in interface ! in :core - module damask_interface ! in :damask_interface:DAMASK_python_interface.f90 + module damask_interface ! in :damask_interface:DAMASK_spectral_interface.f90 - subroutine damask_interface_init(loadcaseParameterIn,geometryParameterIn) ! in :damask_interface:DAMASK_python_interface.f90 + subroutine damask_interface_init(loadcaseParameterIn,geometryParameterIn) ! in :damask_interface:DAMASK_spectral_interface.f90 character(len=1024), intent(in) :: loadcaseParameterIn character(len=1024), intent(in) :: geometryParameterIn end subroutine damask_interface_init diff --git a/code/mesh.f90 b/code/mesh.f90 index a28e257a9..781bc994e 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -261,7 +261,10 @@ mesh_build_subNodeCoords, & mesh_build_ipVolumes, & mesh_build_ipCoordinates, & - mesh_regrid + mesh_regrid, & + mesh_spectral_getResolution, & + mesh_spectral_getDimension, & + mesh_spectral_getHomogenization private :: FE_mapElemtype, & mesh_faceMatch, & mesh_build_FEdata, & @@ -1541,57 +1544,13 @@ enddo ! mesh_Nelems, mesh_Nnodes !******************************************************************** subroutine mesh_spectral_count_nodesAndElements(myUnit) - - use IO, only: IO_lc, & - IO_intValue, & - IO_stringValue, & - IO_stringPos, & - IO_error implicit none integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 7_pInt - integer(pInt), dimension (1+2*maxNchunks) :: myPos - integer(pInt) :: a = 0_pInt, & - b = 0_pInt, & - c = 0_pInt, & - headerLength = 0_pInt, & - i,j - character(len=1024) line,keyword - - mesh_Nnodes = 0_pInt - mesh_Nelems = 0_pInt - - rewind(myUnit) - read(myUnit,'(a1024)') line - myPos = IO_stringPos(line,2_pInt) - keyword = IO_lc(IO_StringValue(line,myPos,2_pInt)) - if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt - else - call IO_error(error_ID=842_pInt) - endif - - rewind(myUnit) - do i = 1_pInt, headerLength - read(myUnit,'(a1024)') line - myPos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_StringValue(line,myPos,1_pInt)) == 'resolution') then - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,myPos,j))) - case('a') - a = IO_intValue(line,myPos,j+1_pInt) - case('b') - b = IO_intValue(line,myPos,j+1_pInt) - case('c') - c = IO_intValue(line,myPos,j+1_pInt) - end select - enddo - mesh_Nelems = a * b * c - mesh_Nnodes = (1_pInt + a)*(1_pInt + b)*(1_pInt + c) - endif - enddo + integer(pInt), dimension(3) :: res + res = mesh_spectral_getResolution(myUnit) + mesh_Nelems = res(1)*res(2)*res(3) + mesh_Nnodes = (1_pInt + res(1))*(1_pInt + res(2))*(1_pInt + res(3)) end subroutine mesh_spectral_count_nodesAndElements @@ -1706,7 +1665,7 @@ end subroutine mesh_abaqus_count_nodesAndElements use IO, only: IO_lc, & IO_stringValue, & IO_stringPos, & - IO_countContinousIntValues + IO_countContinuousIntValues implicit none integer(pInt), intent(in) :: myUnit @@ -1729,7 +1688,7 @@ end subroutine mesh_abaqus_count_nodesAndElements IO_lc(IO_StringValue(line,myPos,2_pInt)) == 'element' ) then mesh_NelemSets = mesh_NelemSets + 1_pInt mesh_maxNelemInSet = max(mesh_maxNelemInSet, & - IO_countContinousIntValues(myUnit)) + IO_countContinuousIntValues(myUnit)) endif enddo @@ -1848,7 +1807,7 @@ subroutine mesh_marc_count_cpElements(myUnit) use IO, only: IO_lc, & IO_stringValue, & IO_stringPos, & - IO_countContinousIntValues + IO_countContinuousIntValues implicit none integer(pInt), intent(in) :: myUnit @@ -1871,7 +1830,7 @@ subroutine mesh_marc_count_cpElements(myUnit) do i=1_pInt,3_pInt+hypoelasticTableStyle ! Skip 3 or 4 lines read (myUnit,610,END=620) line enddo - mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(myUnit) + mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(myUnit) exit endif enddo @@ -1945,7 +1904,7 @@ subroutine mesh_marc_map_elementSets(myUnit) use IO, only: IO_lc, & IO_stringValue, & IO_stringPos, & - IO_continousIntValues + IO_continuousIntValues implicit none integer(pInt), intent(in) :: myUnit @@ -1968,7 +1927,7 @@ subroutine mesh_marc_map_elementSets(myUnit) (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'element' ) ) then elemSet = elemSet+1_pInt mesh_nameElemSet(elemSet) = trim(IO_stringValue(line,myPos,4_pInt)) - mesh_mapElemSet(:,elemSet) = IO_continousIntValues(myUnit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + mesh_mapElemSet(:,elemSet) = IO_continuousIntValues(myUnit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) endif enddo @@ -1986,7 +1945,7 @@ subroutine mesh_abaqus_map_elementSets(myUnit) IO_stringValue, & IO_stringPos, & IO_extractValue, & - IO_continousIntValues, & + IO_continuousIntValues, & IO_error implicit none @@ -2015,7 +1974,7 @@ subroutine mesh_abaqus_map_elementSets(myUnit) if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*elset' ) then elemSet = elemSet + 1_pInt mesh_nameElemSet(elemSet) = trim(IO_extractValue(IO_lc(IO_stringValue(line,myPos,2_pInt)),'elset')) - mesh_mapElemSet(:,elemSet) = IO_continousIntValues(myUnit,mesh_Nelems,mesh_nameElemSet,& + mesh_mapElemSet(:,elemSet) = IO_continuousIntValues(myUnit,mesh_Nelems,mesh_nameElemSet,& mesh_mapElemSet,elemSet-1_pInt) endif enddo @@ -2263,7 +2222,7 @@ subroutine mesh_marc_map_elements(myUnit) use IO, only: IO_lc, & IO_stringValue, & IO_stringPos, & - IO_continousIntValues + IO_continuousIntValues implicit none integer(pInt), intent(in) :: myUnit @@ -2287,7 +2246,7 @@ subroutine mesh_marc_map_elements(myUnit) do i=1_pInt,3_pInt+hypoelasticTableStyle ! skip three (or four if new table style!) lines read (myUnit,610,END=660) line enddo - contInts = IO_continousIntValues(myUnit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) + contInts = IO_continuousIntValues(myUnit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) do i = 1_pInt,contInts(1) cpElem = cpElem+1_pInt mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i) @@ -2517,83 +2476,33 @@ subroutine mesh_abaqus_count_cpSizes(myUnit) !******************************************************************** subroutine mesh_spectral_build_nodes(myUnit) - use IO, only: IO_lc, & - IO_stringValue, & - IO_stringPos, & - IO_error, & - IO_floatValue, & - IO_intValue + use IO, only: & + IO_error implicit none integer(pInt), intent(in) :: myUnit - - integer(pInt), parameter :: maxNchunks = 7_pInt - integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - integer(pInt) :: a = 1_pInt, & - b = 1_pInt, & - c = 1_pInt, & - headerLength = 0_pInt,i,j,n - real(pReal) :: x = 1.0_pReal, & - y = 1.0_pReal, & - z = 1.0_pReal - logical :: gotResolution = .false. ,gotDimension = .false. - character(len=1024) :: line, keyword + integer(pInt) :: n + integer(pInt), dimension(3) :: res = 1_pInt + real(pReal), dimension(3) :: geomdim = 1.0_pReal allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal - - rewind(myUnit) - read(myUnit,'(a1024)') line - myPos = IO_stringPos(line,2_pInt) - keyword = IO_lc(IO_StringValue(line,myPos,2_pInt)) - if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt - else - call IO_error(error_ID=842_pInt) - endif - - rewind(myUnit) - do i = 1_pInt, headerLength - read(myUnit,'(a1024)') line - myPos = IO_stringPos(line,maxNchunks) - select case ( IO_lc(IO_StringValue(line,myPos,1_pInt)) ) - case ('dimension') - gotDimension = .true. - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,myPos,j))) - case('x') - x = IO_floatValue(line,myPos,j+1_pInt) - case('y') - y = IO_floatValue(line,myPos,j+1_pInt) - case('z') - z = IO_floatValue(line,myPos,j+1_pInt) - end select - enddo - case ('resolution') - gotResolution = .true. - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,myPos,j))) - case('a') - a = 1_pInt + IO_intValue(line,myPos,j+1_pInt) - case('b') - b = 1_pInt + IO_intValue(line,myPos,j+1_pInt) - case('c') - c = 1_pInt + IO_intValue(line,myPos,j+1_pInt) - end select - enddo - end select - enddo -! --- sanity checks --- - - if ((.not. gotDimension) .or. (.not. gotResolution)) call IO_error(error_ID=842_pInt) - if ((a < 1_pInt) .or. (b < 1_pInt) .or. (c < 0_pInt)) call IO_error(error_ID=843_pInt) ! 1_pInt is already added - if ((x <= 0.0_pReal) .or. (y <= 0.0_pReal) .or. (z <= 0.0_pReal)) call IO_error(error_ID=844_pInt) + res = mesh_spectral_getResolution(myUnit) + 1_pInt + geomdim = mesh_spectral_getDimension(myUnit) + + if ((res(1) < 1_pInt) .or. (res(2) < 1_pInt) .or. (res(3) < 0_pInt)) & + call IO_error(error_ID=843_pInt) ! 1_pInt is already added + if ((geomdim(1) <= 0.0_pReal) .or. (geomdim(2) <= 0.0_pReal) .or. (geomdim(3) <= 0.0_pReal)) & + call IO_error(error_ID=844_pInt) forall (n = 0_pInt:mesh_Nnodes-1_pInt) - mesh_node0(1,n+1_pInt) = x * real(mod(n,a),pReal) / real(a-1_pInt,pReal) - mesh_node0(2,n+1_pInt) = y * real(mod(n/a,b),pReal) / real(b-1_pInt,pReal) - mesh_node0(3,n+1_pInt) = z * real(mod(n/a/b,c),pReal) / real(c-1_pInt,pReal) + mesh_node0(1,n+1_pInt) = geomdim(1) * real(mod(n,res(1) ),pReal) & + / real(res(1) -1_pInt,pReal) + mesh_node0(2,n+1_pInt) = geomdim(2) * real(mod(n/res(1) ,res(2)),pReal) & + / real(res(2) -1_pInt,pReal) + mesh_node0(3,n+1_pInt) = geomdim(3) * real(mod(n/res(1) /res(2),res(3)),pReal) & + / real(res(3) -1_pInt,pReal) end forall mesh_node = mesh_node0 !why? @@ -2724,25 +2633,27 @@ subroutine mesh_spectral_build_elements(myUnit) use IO, only: IO_lc, & IO_stringValue, & - IO_floatValue, & IO_stringPos, & IO_error, & - IO_continousIntValues, & + IO_continuousIntValues, & IO_intValue, & - IO_countContinousIntValues + IO_countContinuousIntValues implicit none integer(pInt), intent(in) :: myUnit integer(pInt), parameter :: maxNchunks = 7_pInt integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - integer(pInt) :: a = 1_pInt, b = 1_pInt, c = 1_pInt + integer(pInt), dimension(3) :: res integer(pInt) :: e, i, j, homog = 0_pInt, headerLength = 0_pInt, maxIntCount integer(pInt), dimension(:), allocatable :: microstructures integer(pInt), dimension(1,1) :: dummySet = 0_pInt character(len=65536) :: line,keyword character(len=64), dimension(1) :: dummyName = '' + res = mesh_spectral_getResolution(myUnit) + homog = mesh_spectral_getHomogenization(myUnit) + rewind(myUnit) read(myUnit,'(a65536)') line myPos = IO_stringPos(line,2_pInt) @@ -2756,29 +2667,13 @@ subroutine mesh_spectral_build_elements(myUnit) rewind(myUnit) do i = 1_pInt, headerLength read(myUnit,'(a65536)') line - myPos = IO_stringPos(line,maxNchunks) - select case ( IO_lc(IO_StringValue(line,myPos,1_pInt)) ) - case ('resolution') - do j = 2_pInt,6_pInt,2_pInt - select case (IO_lc(IO_stringValue(line,myPos,j))) - case('a') - a = 1_pInt + IO_intValue(line,myPos,j+1_pInt) - case('b') - b = 1_pInt + IO_intValue(line,myPos,j+1_pInt) - case('c') - c = 1_pInt + IO_intValue(line,myPos,j+1_pInt) - end select - enddo - case ('homogenization') - homog = IO_intValue(line,myPos,2_pInt) - end select enddo maxIntCount = 0_pInt i = 1_pInt do while (i > 0_pInt) - i = IO_countContinousIntValues(myUnit) + i = IO_countContinuousIntValues(myUnit) maxIntCount = max(maxIntCount, i) enddo @@ -2792,21 +2687,22 @@ subroutine mesh_spectral_build_elements(myUnit) e = 0_pInt do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) - microstructures = IO_continousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements + microstructures = IO_continuousIntValues(myUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements do i = 1_pInt,microstructures(1_pInt) e = e+1_pInt ! valid element entry mesh_element( 1,e) = e ! FE id mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type mesh_element( 3,e) = homog ! homogenization mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure - mesh_element( 5,e) = e + (e-1_pInt)/(a-1_pInt) + ((e-1_pInt)/((a-1_pInt)*(b-1_pInt)))*a ! base node + mesh_element( 5,e) = e + (e-1_pInt)/(res(1)-1_pInt) + & + ((e-1_pInt)/((res(1)-1_pInt)*(res(2)-1_pInt)))*res(1) ! base node mesh_element( 6,e) = mesh_element(5,e) + 1_pInt - mesh_element( 7,e) = mesh_element(5,e) + a + 1_pInt - mesh_element( 8,e) = mesh_element(5,e) + a - mesh_element( 9,e) = mesh_element(5,e) + a * b ! second floor base node + mesh_element( 7,e) = mesh_element(5,e) + res(1) + 1_pInt + mesh_element( 8,e) = mesh_element(5,e) + res(1) + mesh_element( 9,e) = mesh_element(5,e) + res(1) * res(2) ! second floor base node mesh_element(10,e) = mesh_element(9,e) + 1_pInt - mesh_element(11,e) = mesh_element(9,e) + a + 1_pInt - mesh_element(12,e) = mesh_element(9,e) + a + mesh_element(11,e) = mesh_element(9,e) + res(1) + 1_pInt + mesh_element(12,e) = mesh_element(9,e) + res(1) mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) enddo @@ -2832,7 +2728,7 @@ subroutine mesh_marc_build_elements(myUnit) IO_skipChunks, & IO_stringPos, & IO_intValue, & - IO_continousIntValues + IO_continuousIntValues implicit none integer(pInt), intent(in) :: myUnit @@ -2890,7 +2786,7 @@ subroutine mesh_marc_build_elements(myUnit) read (myUnit,610,END=630) line ! read extra line read (myUnit,610,END=630) line ! read extra line endif - contInts = IO_continousIntValues& ! get affected elements + contInts = IO_continuousIntValues& ! get affected elements (myUnit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) do i = 1_pInt,contInts(1) e = mesh_FEasCP('elem',contInts(1_pInt+i)) @@ -3629,9 +3525,13 @@ deallocate(mesh_HomogMicro) end subroutine mesh_tell_statistics subroutine mesh_regrid(res,resNew) !use new_res=0.0 for automatic determination of new grid - use prec, only: pInt, pReal - use DAMASK_interface, only : getSolverJobName - use IO, only : IO_read_jobBinaryFile + use prec, only: & + pInt, & + pReal + use DAMASK_interface, only: & + getSolverJobName + use IO, only: & + IO_read_jobBinaryFile integer(pInt), dimension(3), intent(in) :: res integer(pInt), dimension(3), intent(inout) :: resNew @@ -3652,5 +3552,186 @@ subroutine mesh_regrid(res,resNew) !use new_res=0.0 for automatic dete end subroutine mesh_regrid +function mesh_spectral_getDimension(myUnit) + use IO, only: & + IO_open_file, & + IO_stringPos, & + IO_lc, & + IO_stringValue, & + IO_intValue, & + IO_floatValue, & + IO_error + use DAMASK_interface, only: & + getModelName, & + InputFileExtension + + implicit none + integer(pInt), dimension(1_pInt + 4_pInt*2_pInt) :: positions ! for a,b c + integer(pInt), optional :: myUnit + integer(pInt) :: headerLength = 0_pInt + real(pReal), dimension(3) :: mesh_spectral_getDimension + character(len=1024) :: line, & + keyword + integer(pInt) :: i, j + logical :: gotDimension = .false. + + if ( .not. present(myUnit)) myUnit = 869 + + call IO_open_file(myUnit,trim(getModelName())//InputFileExtension) + rewind(myUnit) + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,2_pInt) + keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) + if (keyword(1:4) == 'head') then + headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt + else + call IO_error(error_ID=842_pInt) + endif + rewind(myUnit) + do i = 1_pInt, headerLength + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,7_pInt) + select case ( IO_lc(IO_StringValue(line,positions,1)) ) + case ('dimension') + gotDimension = .true. + do j = 2_pInt,6_pInt,2_pInt + select case (IO_lc(IO_stringValue(line,positions,j))) + case('x') + mesh_spectral_getDimension(1) = IO_floatValue(line,positions,j+1_pInt) + case('y') + mesh_spectral_getDimension(2) = IO_floatValue(line,positions,j+1_pInt) + case('z') + mesh_spectral_getDimension(3) = IO_floatValue(line,positions,j+1_pInt) + end select + enddo + end select + enddo + close(myUnit) + if (.not. gotDimension) & + call IO_error(error_ID = 845_pInt, ext_msg='dimension') + if (any(mesh_spectral_getDimension<=0.0_pReal)) & + call IO_error(error_ID = 802_pInt, ext_msg='dimension') + +end function mesh_spectral_getDimension + + +function mesh_spectral_getResolution(myUnit) + use IO, only: & + IO_open_file, & + IO_stringPos, & + IO_lc, & + IO_stringValue, & + IO_intValue, & + IO_floatValue, & + IO_error + use DAMASK_interface, only: & + getModelName, & + InputFileExtension + + implicit none + integer(pInt), dimension(1_pInt + 4_pInt*2_pInt) :: positions ! for a,b c + integer(pInt), optional :: myUnit + integer(pInt) :: headerLength = 0_pInt + integer(pInt), dimension(3) :: mesh_spectral_getResolution + character(len=1024) :: line, & + keyword + integer(pInt) :: i, j + logical :: gotResolution = .false. + + if ( .not. present(myUnit)) myUnit = 869 + + call IO_open_file(myUnit,trim(getModelName())//InputFileExtension) + rewind(myUnit) + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,2_pInt) + keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) + if (keyword(1:4) == 'head') then + headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt + else + call IO_error(error_ID=842_pInt) + endif + rewind(myUnit) + do i = 1_pInt, headerLength + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,7_pInt) + select case ( IO_lc(IO_StringValue(line,positions,1)) ) + case ('resolution') + gotResolution = .true. + do j = 2_pInt,6_pInt,2_pInt + select case (IO_lc(IO_stringValue(line,positions,j))) + case('a') + mesh_spectral_getResolution(1) = IO_intValue(line,positions,j+1_pInt) + case('b') + mesh_spectral_getResolution(2) = IO_intValue(line,positions,j+1_pInt) + case('c') + mesh_spectral_getResolution(3) = IO_intValue(line,positions,j+1_pInt) + end select + enddo + end select + enddo + close(myUnit) + + if (.not. gotResolution) & + call IO_error(error_ID = 845_pInt, ext_msg='resolution') + if(mod(mesh_spectral_getResolution(1),2_pInt)/=0_pInt .or.& + mod(mesh_spectral_getResolution(2),2_pInt)/=0_pInt .or.& + (mod(mesh_spectral_getResolution(3),2_pInt)/=0_pInt .and. & + mesh_spectral_getResolution(3)/= 1_pInt))& + call IO_error(error_ID = 802_pInt, ext_msg='resolution') + +end function mesh_spectral_getResolution + +function mesh_spectral_getHomogenization(myUnit) + use IO, only: & + IO_open_file, & + IO_stringPos, & + IO_lc, & + IO_stringValue, & + IO_intValue, & + IO_error + use DAMASK_interface, only: & + getModelName, & + InputFileExtension + + implicit none + integer(pInt), dimension(1_pInt + 2_pInt*2_pInt) :: positions + integer(pInt), optional :: myUnit + integer(pInt) :: headerLength = 0_pInt + integer(pInt) :: mesh_spectral_getHomogenization + character(len=1024) :: line, & + keyword + integer(pInt) :: i, j + logical :: gotHomogenization = .false. + + if ( .not. present(myUnit)) myUnit = 869 + + call IO_open_file(myUnit,trim(getModelName())//InputFileExtension) + rewind(myUnit) + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,2_pInt) + keyword = IO_lc(IO_StringValue(line,positions,2_pInt)) + if (keyword(1:4) == 'head') then + headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt + else + call IO_error(error_ID=842_pInt) + endif + rewind(myUnit) + do i = 1_pInt, headerLength + read(myUnit,'(a1024)') line + positions = IO_stringPos(line,7_pInt) + select case ( IO_lc(IO_StringValue(line,positions,1)) ) + case ('homogenization') + gotHomogenization = .true. + mesh_spectral_getHomogenization = IO_intValue(line,positions,2_pInt) + end select + enddo + close(myUnit) + + if (.not. gotHomogenization ) & + call IO_error(error_ID = 845_pInt, ext_msg='homogenization') + if (mesh_spectral_getHomogenization<1_pInt) & + call IO_error(error_ID = 802_pInt, ext_msg='homogenization') + +end function mesh_spectral_getHomogenization end module mesh diff --git a/code/numerics.f90 b/code/numerics.f90 index 7b2cac5bd..39c0d1bc4 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -411,7 +411,8 @@ subroutine numerics_init if (err_stress_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolabs') if (itmax <= 1.0_pInt) call IO_error(301_pInt,ext_msg='itmax') if (itmin > itmax) call IO_error(301_pInt,ext_msg='itmin') - + if (update_gamma .and. & + .not. memory_efficient) call IO_error(error_ID = 847_pInt) if (fixedSeed <= 0_pInt) then !$OMP CRITICAL (write2out) write(6,'(a)') ' Random is random!' diff --git a/code/DAMASK_quit.f90 b/code/spectral_quit.f90 similarity index 98% rename from code/DAMASK_quit.f90 rename to code/spectral_quit.f90 index 233ad5b6b..2ed5740b2 100644 --- a/code/DAMASK_quit.f90 +++ b/code/spectral_quit.f90 @@ -46,7 +46,7 @@ subroutine quit(stop_id) integer, dimension(8) :: dateAndTime ! type default integer call date_and_time(values = dateAndTime) - write(6,'(/,a)') 'Terminated on:' + write(6,'(/,a)') 'DAMASK terminated on:' write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',& dateAndTime(2),'/',& dateAndTime(1)