From c786336af374551df180987e23c6a352123b90d2 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Mon, 13 Feb 2012 17:41:27 +0000 Subject: [PATCH] reordered (and partly redistributed) error message identifiers, deleted those which are not in use anymore. all constitutive as well numerics now raises an error, if an unknown keyword is found in the respective config file --- code/CPFEM.f90 | 193 +++++------ code/DAMASK_marc.f90 | 2 +- code/DAMASK_spectral.f90 | 57 ++- code/FEsolving.f90 | 103 +++--- code/IO.f90 | 518 +++++++++++++++------------- code/constitutive.f90 | 10 +- code/constitutive_dislotwin.f90 | 39 ++- code/constitutive_j2.f90 | 89 +++-- code/constitutive_nonlocal.f90 | 138 ++++---- code/constitutive_phenopowerlaw.f90 | 22 +- code/constitutive_titanmod.f90 | 53 +-- code/crystallite.f90 | 6 +- code/debug.f90 | 4 +- code/homogenization.f90 | 12 +- code/lattice.f90 | 8 +- code/material.f90 | 28 +- code/math.f90 | 32 +- code/mesh.f90 | 121 ++++--- code/numerics.f90 | 124 ++++--- 19 files changed, 796 insertions(+), 763 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index e9beec4bb..9c609da6b 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -149,54 +149,54 @@ subroutine CPFEM_init() write(6,'(a)') '<< CPFEM >> Restored state variables of last converged step from binary files' !$OMP END CRITICAL (write2out) endif - if (IO_read_jobBinaryFile(777,'recordedPhase',FEmodelGeometry,size(material_phase))) then - read (777,rec=1) material_phase - close (777) - endif - if (IO_read_jobBinaryFile(777,'convergedF',FEmodelGeometry,size(crystallite_F0))) then - read (777,rec=1) crystallite_F0 - close (777) - endif - if (IO_read_jobBinaryFile(777,'convergedFp',FEmodelGeometry,size(crystallite_Fp0))) then - read (777,rec=1) crystallite_Fp0 - close (777) - endif - if (IO_read_jobBinaryFile(777,'convergedLp',FEmodelGeometry,size(crystallite_Lp0))) then - read (777,rec=1) crystallite_Lp0 - close (777) - endif - if (IO_read_jobBinaryFile(777,'convergeddPdF',FEmodelGeometry,size(crystallite_dPdF0))) then - read (777,rec=1) crystallite_dPdF0 - close (777) - endif - if (IO_read_jobBinaryFile(777,'convergedTstar',FEmodelGeometry,size(crystallite_Tstar0_v))) then - read (777,rec=1) crystallite_Tstar0_v - close (777) - endif - if (IO_read_jobBinaryFile(777,'convergedStateConst',FEmodelGeometry)) then - m = 0_pInt - do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems - do l = 1,size(constitutive_state0(i,j,k)%p) - m = m+1_pInt - read(777,rec=m) constitutive_state0(i,j,k)%p(l) - enddo - enddo; enddo; enddo - close (777) - endif - if (IO_read_jobBinaryFile(777,'convergedStateHomog',FEmodelGeometry)) then - m = 0_pInt - do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips - do l = 1,homogenization_sizeState(j,k) - m = m+1_pInt - read(777,rec=m) homogenization_state0(j,k)%p(l) - enddo - enddo; enddo - close (777) - endif - if (IO_read_jobBinaryFile(777,'convergeddcsdE',FEmodelGeometry,size(CPFEM_dcsdE))) then - read (777,rec=1) CPFEM_dcsdE - close (777) - endif + + call IO_read_jobBinaryFile(777,'recordedPhase',FEmodelGeometry,size(material_phase)) + read (777,rec=1) material_phase + close (777) + + call IO_read_jobBinaryFile(777,'convergedF',FEmodelGeometry,size(crystallite_F0)) + read (777,rec=1) crystallite_F0 + close (777) + + call IO_read_jobBinaryFile(777,'convergedFp',FEmodelGeometry,size(crystallite_Fp0)) + read (777,rec=1) crystallite_Fp0 + close (777) + + call IO_read_jobBinaryFile(777,'convergedLp',FEmodelGeometry,size(crystallite_Lp0)) + read (777,rec=1) crystallite_Lp0 + close (777) + + call IO_read_jobBinaryFile(777,'convergeddPdF',FEmodelGeometry,size(crystallite_dPdF0)) + read (777,rec=1) crystallite_dPdF0 + close (777) + + call IO_read_jobBinaryFile(777,'convergedTstar',FEmodelGeometry,size(crystallite_Tstar0_v)) + read (777,rec=1) crystallite_Tstar0_v + close (777) + + call IO_read_jobBinaryFile(777,'convergedStateConst',FEmodelGeometry) + m = 0_pInt + do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems + do l = 1,size(constitutive_state0(i,j,k)%p) + m = m+1_pInt + read(777,rec=m) constitutive_state0(i,j,k)%p(l) + enddo + enddo; enddo; enddo + close (777) + + call IO_read_jobBinaryFile(777,'convergedStateHomog',FEmodelGeometry) + m = 0_pInt + do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips + do l = 1,homogenization_sizeState(j,k) + m = m+1_pInt + read(777,rec=m) homogenization_state0(j,k)%p(l) + enddo + enddo; enddo + close (777) + + call IO_read_jobBinaryFile(777,'convergeddcsdE',FEmodelGeometry,size(CPFEM_dcsdE)) + read (777,rec=1) CPFEM_dcsdE + close (777) restartRead = .false. endif ! *** end of restoring @@ -427,54 +427,55 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, write(6,'(a)') '<< CPFEM >> Writing state variables of last converged step to binary files' !$OMP END CRITICAL (write2out) endif - if (IO_write_jobBinaryFile(777,'recordedPhase',size(material_phase))) then - write (777,rec=1) material_phase - close (777) - endif - if (IO_write_jobBinaryFile(777,'convergedF',size(crystallite_F0))) then - write (777,rec=1) crystallite_F0 - close (777) - endif - if (IO_write_jobBinaryFile(777,'convergedFp',size(crystallite_Fp0))) then - write (777,rec=1) crystallite_Fp0 - close (777) - endif - if (IO_write_jobBinaryFile(777,'convergedLp',size(crystallite_Lp0))) then - write (777,rec=1) crystallite_Lp0 - close (777) - endif - if (IO_write_jobBinaryFile(777,'convergeddPdF',size(crystallite_dPdF0))) then - write (777,rec=1) crystallite_dPdF0 - close (777) - endif - if (IO_write_jobBinaryFile(777,'convergedTstar',size(crystallite_Tstar0_v))) then - write (777,rec=1) crystallite_Tstar0_v - close (777) - endif - if (IO_write_jobBinaryFile(777,'convergedStateConst')) then - m = 0_pInt - do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems - do l = 1,size(constitutive_state0(i,j,k)%p) - m = m+1_pInt - write(777,rec=m) constitutive_state0(i,j,k)%p(l) - enddo - enddo; enddo; enddo - close (777) - endif - if (IO_write_jobBinaryFile(777,'convergedStateHomog')) then - m = 0_pInt - do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips - do l = 1,homogenization_sizeState(j,k) - m = m+1_pInt - write(777,rec=m) homogenization_state0(j,k)%p(l) - enddo - enddo; enddo - close (777) - endif - if (IO_write_jobBinaryFile(777,'convergeddcsdE',size(CPFEM_dcsdE))) then - write (777,rec=1) CPFEM_dcsdE - close (777) - endif + + call IO_write_jobBinaryFile(777,'recordedPhase',size(material_phase)) + write (777,rec=1) material_phase + close (777) + + call IO_write_jobBinaryFile(777,'convergedF',size(crystallite_F0)) + write (777,rec=1) crystallite_F0 + close (777) + + call IO_write_jobBinaryFile(777,'convergedFp',size(crystallite_Fp0)) + write (777,rec=1) crystallite_Fp0 + close (777) + + call IO_write_jobBinaryFile(777,'convergedLp',size(crystallite_Lp0)) + write (777,rec=1) crystallite_Lp0 + close (777) + + call IO_write_jobBinaryFile(777,'convergeddPdF',size(crystallite_dPdF0)) + write (777,rec=1) crystallite_dPdF0 + close (777) + + call IO_write_jobBinaryFile(777,'convergedTstar',size(crystallite_Tstar0_v)) + write (777,rec=1) crystallite_Tstar0_v + close (777) + + call IO_write_jobBinaryFile(777,'convergedStateConst') + m = 0_pInt + do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems + do l = 1,size(constitutive_state0(i,j,k)%p) + m = m+1_pInt + write(777,rec=m) constitutive_state0(i,j,k)%p(l) + enddo + enddo; enddo; enddo + close (777) + + call IO_write_jobBinaryFile(777,'convergedStateHomog') + m = 0_pInt + do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips + do l = 1,homogenization_sizeState(j,k) + m = m+1_pInt + write(777,rec=m) homogenization_state0(j,k)%p(l) + enddo + enddo; enddo + close (777) + + call IO_write_jobBinaryFile(777,'convergeddcsdE',size(CPFEM_dcsdE)) + write (777,rec=1) CPFEM_dcsdE + close (777) + endif ! * end of dumping endif diff --git a/code/DAMASK_marc.f90 b/code/DAMASK_marc.f90 index 5a76d9df1..0b13b7a87 100644 --- a/code/DAMASK_marc.f90 +++ b/code/DAMASK_marc.f90 @@ -400,7 +400,7 @@ subroutine plotv(& real(pReal) v, t(*) integer(pInt) m, nn, layer, ndi, nshear, jpltcd - if (jpltcd > materialpoint_sizeResults) call IO_error(667,jpltcd) ! complain about out of bounds error + if (jpltcd > materialpoint_sizeResults) call IO_error(700_pInt,jpltcd) ! complain about out of bounds error v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m)) return diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 91eaaf9f9..c36cccf88 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -210,7 +210,7 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! reading the load case file and allocate data structure containing load cases path = getLoadcaseName() - if (.not. IO_open_file(myUnit,path)) call IO_error(error_ID = 30_pInt,ext_msg = trim(path)) + call IO_open_file(myUnit,path) rewind(myUnit) do read(myUnit,'(a1024)',END = 100) line @@ -232,7 +232,7 @@ program DAMASK_spectral 100 N_Loadcases = N_n if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check - call IO_error(error_ID=37_pInt,ext_msg = trim(path)) ! error message for incomplete loadcase + call IO_error(error_ID=837_pInt,ext_msg = trim(path)) ! error message for incomplete loadcase allocate (bc(N_Loadcases)) @@ -310,14 +310,13 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- 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 = 47_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() - if (.not. IO_open_file(myUnit,trim(path)//InputFileExtension))& - call IO_error(error_ID=101_pInt,ext_msg = trim(path)//InputFileExtension) + call IO_open_file(myUnit,trim(path)//InputFileExtension) rewind(myUnit) read(myUnit,'(a1024)') line positions = IO_stringPos(line,2_pInt) @@ -325,7 +324,7 @@ program DAMASK_spectral if (keyword(1:4) == 'head') then headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt else - call IO_error(error_ID=42_pInt) + call IO_error(error_ID=842_pInt) endif rewind(myUnit) @@ -367,12 +366,12 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! sanity checks of geometry parameters if (.not.(gotDimension .and. gotHomogenization .and. gotResolution))& - call IO_error(error_ID = 45_pInt) - if (any(geomdim<=0.0_pReal)) call IO_error(error_ID = 102_pInt) + 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 = 103_pInt) + call IO_error(error_ID = 803_pInt) !-------------------------------------------------------------------------------------------------- ! variables derived from resolution @@ -413,7 +412,7 @@ program DAMASK_spectral 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 = 32_pInt ! each row should be either fully or not at all defined + 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:' else @@ -433,17 +432,17 @@ program DAMASK_spectral print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency - if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 31 ! exclusive or masking only + 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]))) & - errorID = 38_pInt ! no rotation is allowed by stress BC + errorID = 838_pInt ! no rotation is allowed by stress BC if (any(abs(math_mul33x33(bc(loadcase)%rotation,math_transpose33(bc(loadcase)%rotation))& -math_I3) > reshape(spread(rotation_tol,1,9),[ 3,3]))& .or. abs(math_det33(bc(loadcase)%rotation)) > 1.0_pReal + rotation_tol)& - errorID = 46_pInt ! given rotation matrix contains strain - if (bc(loadcase)%time < 0.0_pReal) errorID = 34_pInt ! negative time increment - if (bc(loadcase)%incs < 1_pInt) errorID = 35_pInt ! non-positive incs count - if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 36_pInt ! non-positive result frequency + errorID = 846_pInt ! given rotation matrix contains strain + if (bc(loadcase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment + if (bc(loadcase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count + if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) enddo @@ -483,7 +482,7 @@ program DAMASK_spectral c0_reference = c_current * wgt ! linear reference material stiffness c0_66 = math_Mandel3333to66(c0_reference) call math_invert(6_pInt, c0_66, s0_66, i, errmatinv) ! invert in mandel notation - if(errmatinv) call IO_error(error_ID=800_pInt) + if(errmatinv) call IO_error(error_ID=400_pInt) s0_reference = math_Mandel66to3333(s0_66) !-------------------------------------------------------------------------------------------------- @@ -491,11 +490,10 @@ program DAMASK_spectral if (restartInc > 1_pInt) then ! using old values from file if (debugRestart) print '(a,i6,a)' , 'Reading values of increment ',& restartInc - 1_pInt,' from file' - if (IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& - trim(getSolverJobName()),size(defgrad))) then - read (777,rec=1) defgrad - close (777) - endif + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& + trim(getSolverJobName()),size(defgrad)) + read (777,rec=1) defgrad + close (777) defgradold = defgrad defgradAim = 0.0_pReal do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) @@ -544,11 +542,11 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! 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=108_pInt) ! check for correct precision in C + 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 = 109_pInt) + if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt) call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) endif #endif @@ -723,7 +721,7 @@ program DAMASK_spectral c_reduced(k,j) = c_prev99(n,m) endif; enddo; endif; enddo call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) ! invert reduced stiffness - if(errmatinv) call IO_error(error_ID=800_pInt) + if(errmatinv) call IO_error(error_ID=400_pInt) s_prev99 = 0.0_pReal ! build full compliance k = 0_pInt do n = 1_pInt,9_pInt @@ -1097,10 +1095,9 @@ program DAMASK_spectral mod(inc - 1_pInt,bc(loadcase)%restartFrequency) == 0_pInt) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?) restartWrite = .true. print '(A)', 'writing converged results for restart' - if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then ! writing deformation gradient field to file - write (777,rec=1) defgrad - close (777) - endif + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad)) ! writing deformation gradient field to file + write (777,rec=1) defgrad + close (777) restartInc=totalIncsCounter endif @@ -1111,7 +1108,7 @@ program DAMASK_spectral !c0_99 = math_Plain3333to99(c0_reference) ! call math_invert(9_pInt, s0_99, c0_99, i, errmatinv) ! invert reduced stiffness - ! if(errmatinv) call IO_error(error_ID=800_pInt) + ! if(errmatinv) call IO_error(error_ID=400_pInt) ! print*, (c0_reference - math_Plain99to3333(c0_99))/c0_reference ! pause endif diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index 58205297e..782225e1a 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -60,63 +60,61 @@ integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions character(len=1024) line, commandLine FEmodelGeometry = getModelName() - if (IO_open_inputFile(fileunit,FEmodelGeometry)) then - if (trim(FEsolver) == 'Spectral') then - call get_command(commandLine) ! may contain uppercase - do i=1,len(commandLine) - if(64 < iachar(commandLine(i:i)) .and. iachar(commandLine(i:i)) < 91) & - commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase - enddo - if (index(commandLine,'-r ',.true.)>0) & ! look for -r - start = index(commandLine,'-r ',.true.) + 3_pInt ! set to position after trailing space - if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart - start = index(commandLine,'--restart ',.true.) + 10_pInt ! set to position after trailing space + call IO_open_inputFile(fileunit,FEmodelGeometry) + if (trim(FEsolver) == 'Spectral') then + call get_command(commandLine) ! may contain uppercase + do i=1,len(commandLine) + if(64 < iachar(commandLine(i:i)) .and. iachar(commandLine(i:i)) < 91) & + commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase + enddo + if (index(commandLine,'-r ',.true.)>0) & ! look for -r + start = index(commandLine,'-r ',.true.) + 3_pInt ! set to position after trailing space + if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart + start = index(commandLine,'--restart ',.true.) + 10_pInt ! set to position after trailing space - if(start /= 0_pInt) then ! found something - length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument? - read(commandLine(start:start+length),'(I12)') restartInc ! read argument - restartRead = restartInc > 0_pInt - if(restartInc <= 0_pInt) then - call IO_warning(warning_ID=34_pInt) - restartInc = 1_pInt - endif + if(start /= 0_pInt) then ! found something + length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument? + read(commandLine(start:start+length),'(I12)') restartInc ! read argument + restartRead = restartInc > 0_pInt + if(restartInc <= 0_pInt) then + call IO_warning(warning_ID=34_pInt) + restartInc = 1_pInt endif - else - rewind(fileunit) - do - read (fileunit,'(a1024)',END=100) line - positions = IO_stringPos(line,maxNchunks) - tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key - select case(tag) - case ('solver') - read (fileunit,'(a1024)',END=100) line ! next line - positions = IO_stringPos(line,maxNchunks) - symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt) - case ('restart') - read (fileunit,'(a1024)',END=100) line ! next line - positions = IO_stringPos(line,maxNchunks) - restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt - restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt - case ('*restart') - do i=2,positions(1) - restartWrite = (IO_lc(IO_StringValue(line,positions,i)) == 'write') .or. restartWrite - restartRead = (IO_lc(IO_StringValue(line,positions,i)) == 'read') .or. restartRead - enddo - if(restartWrite) then - do i=2,positions(1) - restartWrite = (IO_lc(IO_StringValue(line,positions,i)) /= 'frequency=0') .and. restartWrite - enddo - endif - end select - enddo endif else - call IO_error(101_pInt, ext_msg=FEmodelGeometry) ! cannot open input file + rewind(fileunit) + do + read (fileunit,'(a1024)',END=100) line + positions = IO_stringPos(line,maxNchunks) + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + select case(tag) + case ('solver') + read (fileunit,'(a1024)',END=100) line ! next line + positions = IO_stringPos(line,maxNchunks) + symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt) + case ('restart') + read (fileunit,'(a1024)',END=100) line ! next line + positions = IO_stringPos(line,maxNchunks) + restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt + restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt + case ('*restart') + do i=2,positions(1) + restartWrite = (IO_lc(IO_StringValue(line,positions,i)) == 'write') .or. restartWrite + restartRead = (IO_lc(IO_StringValue(line,positions,i)) == 'read') .or. restartRead + enddo + if(restartWrite) then + do i=2,positions(1) + restartWrite = (IO_lc(IO_StringValue(line,positions,i)) /= 'frequency=0') .and. restartWrite + enddo + endif + end select + enddo endif 100 close(fileunit) if (restartRead) then - if(FEsolver == 'Marc' .and. IO_open_logFile(fileunit)) then + if(FEsolver == 'Marc') then + call IO_open_logFile(fileunit) rewind(fileunit) do read (fileunit,'(a1024)',END=200) line @@ -127,7 +125,8 @@ integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions IO_lc(IO_stringValue(line,positions,4_pInt)) == 'id' ) & FEmodelGeometry = IO_StringValue(line,positions,6_pInt) enddo - elseif (FEsolver == 'Abaqus' .and. IO_open_inputFile(fileunit,FEmodelGeometry)) then + elseif (FEsolver == 'Abaqus') then + call IO_open_inputFile(fileunit,FEmodelGeometry) rewind(fileunit) do read (fileunit,'(a1024)',END=200) line @@ -138,10 +137,6 @@ integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions FEmodelGeometry = IO_StringValue(line,positions,1_pInt) endif enddo - elseif (FEsolver == 'Spectral') then - !do nothing - else - call IO_error(106_pInt) ! cannot open file for old job info endif endif diff --git a/code/IO.f90 b/code/IO.f90 index f1145f10d..21ff74327 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -25,9 +25,9 @@ CONTAINS !--------------------------- ! function IO_abaqus_assembleInputFile -! function IO_open_file(unit,relPath) -! function IO_open_inputFile(unit, model) -! function IO_open_logFile(unit) +! subroutine IO_open_file(unit,relPath) +! subroutine IO_open_inputFile(unit, model) +! subroutine IO_open_logFile(unit) ! function IO_hybridIA(Nast,ODFfileName) ! private function hybridIA_reps(dV_V,steps,C) ! function IO_stringPos(line,maxN) @@ -153,21 +153,49 @@ end function ! open existing file to given unit ! path to file is relative to working directory !******************************************************************** - logical function IO_open_file(unit,relPath) + logical function IO_open_file_stat(unit,relPath) use prec, only: pInt use DAMASK_interface implicit none - character(len=*) relPath - integer(pInt) unit + integer(pInt), intent(in) :: unit + character(len=*), intent(in) :: relPath + integer stat + character path + + IO_open_file_stat = .false. + path = trim(getSolverWorkingDirectoryName())//relPath + open(unit,status='old',iostat=stat,file=path) + if (stat == 0) then + IO_open_file_stat = .true. + endif + + endfunction - IO_open_file = .false. + +!******************************************************************** +! open existing file to given unit +! path to file is relative to working directory +!******************************************************************** + subroutine IO_open_file(unit,relPath) + + use prec, only: pInt + use DAMASK_interface - open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//relPath) - IO_open_file = .true. + implicit none + integer(pInt), intent(in) :: unit + character(len=*), intent(in) :: relPath + character path + integer(pInt) stat -100 endfunction + path = trim(getSolverWorkingDirectoryName())//relPath + open(unit,status='old',iostat=stat,file=path) + if (stat /= 0) then + call IO_error(100_pInt,ext_msg=path) + endif + + endsubroutine !******************************************************************** @@ -176,59 +204,74 @@ end function ! : changed the function to open *.inp_assembly, which is basically ! the input file without comment lines and possibly assembled includes !******************************************************************** - logical function IO_open_inputFile(unit,model) + subroutine IO_open_inputFile(unit,model) use prec, only: pReal, pInt use DAMASK_interface implicit none integer(pInt), intent(in) :: unit - character(len=*) model + character(len=*), intent(in) :: model + integer(pInt) stat + character path - IO_open_inputFile = .false. if (FEsolver == 'Abaqus') then - open(unit+1,status='old',err=100,& - file=trim(getSolverWorkingDirectoryName())//& - trim(model)//InputFileExtension) - open(unit,err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(model)//InputFileExtension//'_assembly') - IO_open_inputFile = IO_abaqus_assembleInputFile(unit,unit+1_pInt) ! strip comments and concatenate any "include"s + path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension + open(unit+1,status='old',iostat=stat,file=path) + if (stat /= 0) then + call IO_error(100_pInt,ext_msg=path) + endif + + path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension//'_assembly' + open(unit,iostat=stat,file=path) + if (stat /= 0) then + call IO_error(100_pInt,ext_msg=path) + endif + + if (IO_abaqus_assembleInputFile(unit,unit+1_pInt)) then ! strip comments and concatenate any "include"s + call IO_error(103_pInt) + endif close(unit+1_pInt) + else - open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(model)//InputFileExtension) - IO_open_inputFile = .true. + path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension + open(unit,status='old',iostat=stat,file=path) + if (stat /= 0) then + call IO_error(100_pInt,ext_msg=path) + endif endif -100 endfunction + endsubroutine !******************************************************************** ! open FEM logfile to given unit !******************************************************************** - logical function IO_open_logFile(unit) + subroutine IO_open_logFile(unit) use prec, only: pReal, pInt use DAMASK_interface implicit none + character path + integer(pInt) stat integer(pInt), intent(in) :: unit - IO_open_logFile = .false. - - open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(getSolverJobName())//LogFileExtension) - IO_open_logFile = .true. + path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//LogFileExtension + open(unit,status='old',iostat=stat,file=path) + if (stat /= 0) then + call IO_error(100_pInt,ext_msg=path) + endif -100 endfunction + endsubroutine !******************************************************************** ! open (write) file related to current job ! but with different extension to given unit !******************************************************************** - logical function IO_open_jobFile(unit,newExt) + logical function IO_open_jobFile_stat(unit,newExt) use prec, only: pReal, pInt use DAMASK_interface @@ -236,21 +279,24 @@ end function integer(pInt), intent(in) :: unit character(*), intent(in) :: newExt + character path + integer stat - IO_open_jobFile = .false. + IO_open_jobFile_stat = .false. + path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt + open(unit,status='old',iostat=stat,file=path) + if (stat == 0) then + IO_open_jobFile_stat = .true. + endif - open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(getSolverJobName())//'.'//newExt) - IO_open_jobFile = .true. - -100 endfunction + endfunction !******************************************************************** ! open (write) file related to current job ! but with different extension to given unit !******************************************************************** - logical function IO_write_jobFile(unit,newExt) + subroutine IO_open_jobFile(unit,newExt) use prec, only: pReal, pInt use DAMASK_interface @@ -258,21 +304,47 @@ end function integer(pInt), intent(in) :: unit character(*), intent(in) :: newExt + character path + integer(pInt) stat - IO_write_jobFile = .false. + path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt + open(unit,status='old',iostat=stat,file=path) + if (stat /= 0) then + call IO_error(100_pInt,ext_msg=path) + endif - open(unit,status='replace',err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(getSolverJobName())//'.'//newExt) - IO_write_jobFile = .true. + endsubroutine + + +!******************************************************************** +! open (write) file related to current job +! but with different extension to given unit +!******************************************************************** + subroutine IO_write_jobFile(unit,newExt) + + use prec, only: pReal, pInt + use DAMASK_interface + implicit none + + integer(pInt), intent(in) :: unit + character(*), intent(in) :: newExt + character path + integer(pInt) stat + + path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt + open(unit,status='replace',iostat=stat,file=path) + if (stat /= 0) then + call IO_error(100_pInt,ext_msg=path) + endif -100 endfunction + endsubroutine !******************************************************************** ! open (write) binary file related to current job ! but with different extension to given unit !******************************************************************** - logical function IO_write_jobBinaryFile(unit,newExt,recMultiplier) + subroutine IO_write_jobBinaryFile(unit,newExt,recMultiplier) use prec, only: pReal, pInt use DAMASK_interface @@ -281,27 +353,27 @@ end function integer(pInt), intent(in) :: unit integer(pInt), intent(in), optional :: recMultiplier character(*), intent(in) :: newExt + character path + integer(pInt) stat - IO_write_jobBinaryFile = .false. + path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt if (present(recMultiplier)) then - open(unit,status='replace',form='unformatted',access='direct',recl=pReal*recMultiplier, & - err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(getSolverJobName())//'.'//newExt) - else - open(unit,status='replace',form='unformatted',access='direct',recl=pReal, & - err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(getSolverJobName())//'.'//newExt) + open(unit,status='replace',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path) + else + open(unit,status='replace',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path) + endif + if (stat /= 0) then + call IO_error(100_pInt,ext_msg=path) endif - IO_write_jobBinaryFile = .true. -100 endfunction + endsubroutine !******************************************************************** ! open (read) binary file related to restored job ! and with different extension to given unit !******************************************************************** - logical function IO_read_jobBinaryFile(unit,newExt,jobName,recMultiplier) + subroutine IO_read_jobBinaryFile(unit,newExt,jobName,recMultiplier) use prec, only: pReal, pInt use DAMASK_interface @@ -310,20 +382,20 @@ end function integer(pInt), intent(in) :: unit integer(pInt), intent(in), optional :: recMultiplier character(*), intent(in) :: newExt, jobName + character path + integer(pInt) stat - IO_read_jobBinaryFile = .false. + path = trim(getSolverWorkingDirectoryName())//trim(jobName)//'.'//newExt if (present(recMultiplier)) then - open(unit,status='old',form='unformatted',access='direct',recl=pReal*recMultiplier, & - err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(jobName)//'.'//newExt) - else - open(unit,status='old',form='unformatted',access='direct',recl=pReal, & - err=100,file=trim(getSolverWorkingDirectoryName())//& - trim(jobName)//'.'//newExt) + open(unit,status='old',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path) + else + open(unit,status='old',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path) + endif + if (stat /= 0) then + call IO_error(100_pInt,ext_msg=path) endif - IO_read_jobBinaryFile = .true. -100 endfunction + endsubroutine !******************************************************************** @@ -374,7 +446,7 @@ end function real(pReal), dimension(:,:,:), allocatable :: dV_V real(pReal), dimension(3,Nast) :: IO_hybridIA - if (.not. IO_open_file(999_pInt,ODFfileName)) goto 100 + call IO_open_file(999_pInt,ODFfileName) !--- parse header of ODF file --- !--- limits in phi1, Phi, phi2 --- @@ -1135,225 +1207,172 @@ endfunction character(len=1024) msg select case (error_ID) - case (30_pInt) - msg = 'could not open spectral loadcase' - case (31_pInt) - msg = 'mask consistency violated in spectral loadcase' - case (32_pInt) - msg = 'ill-defined L (each line should be either fully or not at all defined) in spectral loadcase' - case (34_pInt) - msg = 'negative time increment in spectral loadcase' - case (35_pInt) - msg = 'non-positive increments in spectral loadcase' - case (36_pInt) - msg = 'non-positive result frequency in spectral loadcase' - case (37_pInt) - msg = 'incomplete loadcase' - case (38_pInt) - msg = 'mixed boundary conditions allow rotation' - case (39_pInt) - msg = 'non-positive restart frequency in spectral loadcase' - case (40_pInt) - msg = 'path rectification error' - case (41_pInt) - msg = 'path too long' - case (42_pInt) - msg = 'missing header info in spectral mesh' - case (43_pInt) - msg = 'resolution in spectral mesh' - case (44_pInt) - msg = 'dimension in spectral mesh' - case (45_pInt) - msg = 'incomplete information in spectral mesh header' - case (46_pInt) - msg = 'not a rotation defined for loadcase rotation' - case (47_pInt) - msg = 'updating of gamma operator not possible if it is pre calculated' - case (50_pInt) - msg = 'writing constitutive output description' + + !* file handling errors + case (100_pInt) - msg = 'opening material configuration' + msg = 'could not open file:' case (101_pInt) - msg = 'opening input file' + msg = 'write error for file:' case (102_pInt) - msg = 'non-positive dimension' + msg = 'could not read file:' case (103_pInt) - msg = 'odd resolution given' - case (105_pInt) - msg = 'reading from ODF file' - case (106_pInt) - msg = 'reading info on old job' - case (107_pInt) - msg = 'writing spectralOut file' - case (108_pInt) - msg = 'precistion not suitable for FFTW' - case (109_pInt) - msg = 'initializing FFTW' - case (110_pInt) - msg = 'no homogenization specified via State Variable 2' - case (120_pInt) - msg = 'no microstructure specified via State Variable 3' - case (125_pInt) - msg = 'no entries in config part' - case (130_pInt) - msg = 'homogenization index out of bounds' - case (140_pInt) - msg = 'microstructure index out of bounds' + msg = 'could not assemble input files' + + + !* material error messages and related messages in mesh + case (150_pInt) msg = 'crystallite index out of bounds' - case (155_pInt) + case (151_pInt) msg = 'phase index out of bounds' - case (160_pInt) + case (152_pInt) msg = 'texture index out of bounds' - case (170_pInt) + case (153_pInt) msg = 'sum of phase fractions differs from 1' + case (154_pInt) + msg = 'homogenization index out of bounds' + case (155_pInt) + msg = 'microstructure index out of bounds' + case (156_pInt) + msg = 'reading from ODF file' + case (160_pInt) + msg = 'no entries in config part' + case (170_pInt) + msg = 'no homogenization specified via State Variable 2' case (180_pInt) - msg = 'mismatch of microstructure count and a*b*c in geom file' + msg = 'no microstructure specified via State Variable 3' + + + !* constitutive error messages + case (200_pInt) msg = 'unknown constitution specified' - case (201_pInt) - msg = 'unknown homogenization specified' case (205_pInt) msg = 'unknown lattice structure encountered' + case (210_pInt) - msg = 'negative initial resistance' + msg = 'unknown material parameter for j2 constitutive phase:' case (211_pInt) - msg = 'non-positive reference shear rate' + msg = 'material parameter for j2 constitutive phase out of bounds:' case (212_pInt) - msg = 'non-positive stress exponent' - case (213_pInt) - msg = 'non-positive saturation stress' - case (214_pInt) - msg = 'zero hardening exponent' + msg = 'unknown constitutive output for j2 constitution:' + case (220_pInt) - msg = 'negative initial dislocation density' + msg = 'unknown material parameter for phenopowerlaw constitutive phase:' case (221_pInt) - msg = 'negative Burgers vector' + msg = 'material parameter for phenopowerlaw constitutive phase out of bounds:' case (222_pInt) - msg = 'negative activation energy for edge dislocation glide' - case (223_pInt) - msg = 'zero stackin fault energy' - case (225_pInt) - msg = 'no slip systems specified' - case (226_pInt) - msg = 'non-positive prefactor for dislocation velocity' - case (227_pInt) - msg = 'non-positive prefactor for mean free path' - case (228_pInt) - msg = 'non-positive minimum stable dipole distance' - case (229_pInt) - msg = 'non-positive hardening interaction coefficients' + msg = 'unknown constitutive output for phenopowerlaw constitution:' + case (230_pInt) - msg = 'non-positive atomic volume' + msg = 'unknown material parameter for titanmod constitutive phase:' case (231_pInt) - msg = 'non-positive prefactor for self-diffusion coefficient' + msg = 'material parameter for titanmod constitutive phase out of bounds:' case (232_pInt) - msg = 'non-positive activation energy for self-diffusion' - case (233_pInt) - msg = 'non-positive relevant dislocation density' - case (234_pInt) - msg = 'error in shear banding input' - case (235_pInt) - msg = 'material parameter for nonlocal constitutive phase out of bounds:' - case (236_pInt) - msg = 'unknown material parameter for nonlocal constitutive phase:' - case (237_pInt) - msg = 'unknown constitutive output for nonlocal constitution:' + msg = 'unknown constitutive output for titanmod constitution:' + case (240_pInt) - msg = 'non-positive Taylor factor' + msg = 'unknown material parameter for dislotwin constitutive phase:' case (241_pInt) - msg = 'non-positive hardening exponent' + msg = 'material parameter for dislotwin constitutive phase out of bounds:' case (242_pInt) - msg = 'non-positive relevant slip resistance' - case (260_pInt) - msg = 'non-positive relevant strain' - case (261_pInt) - msg = 'frequency of stiffness update < 0' - case (262_pInt) - msg = 'frequency of Jacobian update of Lp residuum < 0' - case (263_pInt) - msg = 'non-positive perturbation value' - case (264_pInt) - msg = 'limit for homogenization loop too small' - case (265_pInt) - msg = 'limit for crystallite loop too small' - case (266_pInt) - msg = 'limit for state loop too small' - case (267_pInt) - msg = 'limit for stress loop too small' - case (268_pInt) - msg = 'non-positive minimum substep size' - case (269_pInt) - msg = 'non-positive relative state tolerance' - case (270_pInt) - msg = 'Non-positive relative stress tolerance' - case (271_pInt) - msg = 'Non-positive absolute stress tolerance' + msg = 'unknown constitutive output for dislotwin constitution:' + case (243_pInt) + msg = 'zero stackin fault energy' -!* Error messages related to RGC numerical parameters <<>> - case (272_pInt) - msg = 'non-positive relative tolerance of residual in RGC' - case (273_pInt) - msg = 'non-positive absolute tolerance of residual in RGC' - case (274_pInt) - msg = 'non-positive relative maximum of residual in RGC' - case (275_pInt) - msg = 'non-positive absolute maximum of residual in RGC' - case (276_pInt) - msg = 'non-positive penalty perturbation in RGC' - case (277_pInt) - msg = 'non-positive relevant mismatch in RGC' - case (278_pInt) - msg = 'non-positive definite viscosity model in RGC' - case (288_pInt) - msg = 'non-positive maximum threshold of relaxation change in RGC' - case (289_pInt) - msg = 'non-positive definite volume discrepancy penalty in RGC' + case (250_pInt) + msg = 'unknown material parameter for nonlocal constitutive phase:' + case (251_pInt) + msg = 'material parameter for nonlocal constitutive phase out of bounds:' + case (252_pInt) + msg = 'unknown constitutive output for nonlocal constitution:' + case (253_pInt) + msg = 'element type not supported for nonlocal constitution' - case (294_pInt) - msg = 'non-positive tolerance for deformation gradient' - - case (298_pInt) - msg = 'chosen integration method does not exist' - case (299_pInt) - msg = 'chosen perturbation method does not exist' + + !* numerics error messages case (300_pInt) - msg = 'this material can only be used with elements with three direct stress components' - case (500_pInt) - msg = 'unknown lattice type specified' - case (550_pInt) + msg = 'unknown numerics parameter:' + case (301_pInt) + msg = 'numerics parameter out of bounds:' + + + !* math errors + + case (400_pInt) + msg = 'matrix inversion error' + case (401_pInt) + msg = 'math_check: quat -> axisAngle -> quat failed' + case (402_pInt) + msg = 'math_check: quat -> R -> quat failed' + case (403_pInt) + msg = 'math_check: quat -> euler -> quat failed' + case (404_pInt) + msg = 'math_check: R -> euler -> R failed' + case (405_pInt) + msg = 'I_TO_HALTON-error: An input base BASE is <= 1' + case (406_pInt) + msg = 'Prime-error: N must be between 0 and PRIME_MAX' + case (450_pInt) msg = 'unknown symmetry type specified' - case (600_pInt) - msg = 'convergence not reached' - case (610_pInt) - msg = 'stress loop not converged' - case (666_pInt) - msg = 'memory leak detected' - case (667_pInt) + + !* homogenization errors + + case (500_pInt) + msg = 'unknown homogenization specified' + + + !* DAMASK_marc errors + + case (700_pInt) msg = 'invalid materialpoint result requested' - case (670_pInt) - msg = 'math_check: quat -> axisAngle -> quat failed' - case (671_pInt) - msg = 'math_check: quat -> R -> quat failed' - case (672_pInt) - msg = 'math_check: quat -> euler -> quat failed' - case (673_pInt) - msg = 'math_check: R -> euler -> R failed' - case (700_pInt) - msg = 'singular matrix in stress iteration' + !* errors related to spectral solver - case (800_pInt) - msg = 'matrix inversion error' - case (801_pInt) - msg = 'I_TO_HALTON-error: An input base BASE is <= 1' case (802_pInt) - msg = 'Prime-error: N must be between 0 and PRIME_MAX' + msg = 'non-positive dimension' + case (803_pInt) + msg = 'odd resolution given' + case (808_pInt) + msg = 'precision not suitable for FFTW' + case (809_pInt) + msg = 'initializing FFTW' + case (831_pInt) + msg = 'mask consistency violated in spectral loadcase' + case (832_pInt) + msg = 'ill-defined L (each line should be either fully or not at all defined) in spectral loadcase' + case (834_pInt) + msg = 'negative time increment in spectral loadcase' + case (835_pInt) + msg = 'non-positive increments in spectral loadcase' + case (836_pInt) + msg = 'non-positive result frequency in spectral loadcase' + case (837_pInt) + msg = 'incomplete loadcase' + case (838_pInt) + msg = 'mixed boundary conditions allow rotation' + case (842_pInt) + msg = 'missing header info in spectral mesh' + case (843_pInt) + msg = 'resolution in spectral mesh' + case (844_pInt) + msg = 'dimension in spectral mesh' + case (845_pInt) + msg = 'incomplete information in spectral mesh header' + case (846_pInt) + msg = 'not a rotation defined for loadcase rotation' + case (847_pInt) + msg = 'updating of gamma operator not possible if it is pre calculated' + case (880_pInt) + msg = 'mismatch of microstructure count and a*b*c in geom file' + + + !* Error messages related to parsing of Abaqus input file -! Error messages related to parsing of Abaqus input file case (900_pInt) msg = 'PARSE ERROR: Improper definition of nodes in input file (Nnodes < 2)' case (901_pInt) @@ -1378,8 +1397,13 @@ endfunction msg = 'PARSE ERROR: Incorrect element type mapping in ' + !* general error messages + + case (666_pInt) + msg = 'memory leak detected' case default msg = 'Unknown error number...' + end select !$OMP CRITICAL (write2out) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 4db8d12d0..2d6135b46 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -71,7 +71,7 @@ subroutine constitutive_init() use prec, only: pReal,pInt use debug, only: debug_verbosity, debug_selectiveDebugger, debug_e, debug_i, debug_g use numerics, only: numerics_integrator -use IO, only: IO_error, IO_open_file, IO_open_jobFile, IO_write_jobFile +use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips use material use constitutive_j2 @@ -100,8 +100,8 @@ logical knownConstitution ! --- PARSE CONSTITUTIONS FROM CONFIG FILE --- -if (.not. IO_open_jobFile(fileunit,material_localFileExt)) then ! no local material configuration present... - if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file +if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present... + call IO_open_file(fileunit,material_configFile) ! ... open material.config file endif call constitutive_j2_init(fileunit) call constitutive_phenopowerlaw_init(fileunit) @@ -113,9 +113,7 @@ close(fileunit) ! --- WRITE DESCRIPTION FILE FOR CONSTITUTIVE PHASE OUTPUT --- -if(.not. IO_write_jobFile(fileunit,'outputConstitutive')) then ! problems in writing file - call IO_error (50_pInt) -endif +call IO_write_jobFile(fileunit,'outputConstitutive') do p = 1,material_Nphase i = phase_constitutionInstance(p) ! which instance of a constitution is present phase knownConstitution = .true. ! assume valid diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index f0b651aa8..37de1eb92 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -393,6 +393,8 @@ do ! read thru sections of constitutive_dislotwin_sbResistance(i) = IO_floatValue(line,positions,2_pInt) case ('shearbandvelocity') constitutive_dislotwin_sbVelocity(i) = IO_floatValue(line,positions,2_pInt) + case default + call IO_error(240_pInt,ext_msg=tag) end select endif enddo @@ -403,32 +405,31 @@ enddo myStructure = constitutive_dislotwin_structure(i) !* Sanity checks - if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt) - if (sum(constitutive_dislotwin_Nslip(:,i)) <= 0_pInt) call IO_error(225_pInt) - if (sum(constitutive_dislotwin_Ntwin(:,i)) < 0_pInt) call IO_error(225_pInt) !*** + if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt,e=i) + if (sum(constitutive_dislotwin_Nslip(:,i)) <= 0_pInt) call IO_error(241_pInt,e=i,ext_msg='nslip') + if (sum(constitutive_dislotwin_Ntwin(:,i)) < 0_pInt) call IO_error(241_pInt,e=i,ext_msg='ntwin') do f = 1,lattice_maxNslipFamily if (constitutive_dislotwin_Nslip(f,i) > 0_pInt) then - if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(220_pInt) - if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(220_pInt) - if (constitutive_dislotwin_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(221_pInt) - if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt) + if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='rhoEdge0') + if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='rhoEdgeDip0') + if (constitutive_dislotwin_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='slipburgers') + if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='v0') endif enddo do f = 1,lattice_maxNtwinFamily if (constitutive_dislotwin_Ntwin(f,i) > 0_pInt) then - if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221_pInt) !*** - if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(226_pInt) !*** + if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='twinburgers') + if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='ndot0') endif enddo -! if (any(constitutive_dislotwin_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) call IO_error(229) - if (constitutive_dislotwin_CAtomicVolume(i) <= 0.0_pReal) call IO_error(230_pInt) - if (constitutive_dislotwin_D0(i) <= 0.0_pReal) call IO_error(231_pInt) - if (constitutive_dislotwin_Qsd(i) <= 0.0_pReal) call IO_error(232_pInt) - if (constitutive_dislotwin_aTolRho(i) <= 0.0_pReal) call IO_error(233_pInt) - if (constitutive_dislotwin_sbResistance(i) <= 0.0_pReal) call IO_error(234_pInt) - if (constitutive_dislotwin_sbVelocity(i) < 0.0_pReal) call IO_error(235_pInt) - if (constitutive_dislotwin_SFE_0K(i) == 0.0_pReal .and. & - constitutive_dislotwin_dSFE_dT(i) == 0.0_pReal) call IO_error(223_pInt) + if (constitutive_dislotwin_CAtomicVolume(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='cAtomicVolume') + if (constitutive_dislotwin_D0(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='D0') + if (constitutive_dislotwin_Qsd(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='Qsd') + if (constitutive_dislotwin_aTolRho(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='aTolRho') + if (constitutive_dislotwin_sbResistance(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='sbResistance') + if (constitutive_dislotwin_sbVelocity(i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='sbVelocity') + if (constitutive_dislotwin_SFE_0K(i) == 0.0_pReal .AND. & + constitutive_dislotwin_dSFE_dT(i) == 0.0_pReal) call IO_error(243_pInt,e=i) !* Determine total number of active slip or twin systems constitutive_dislotwin_Nslip(:,i) = min(lattice_NslipSystem(:,myStructure),constitutive_dislotwin_Nslip(:,i)) @@ -534,7 +535,7 @@ do i = 1_pInt,maxNinstance case('sb_eigenvectors') mySize = 9_pInt case default - mySize = 0_pInt + call IO_error(242_pInt,ext_msg=constitutive_dislotwin_output(o,i)) end select if (mySize > 0_pInt) then ! any meaningful output found diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 127451dab..4eb8c5d56 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -83,15 +83,14 @@ subroutine constitutive_j2_init(file) !************************************** !* Module initialization * !************************************** - 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, pReal use math, only: math_Mandel3333to66, math_Voigt66to3333 use IO use material use debug, only: debug_verbosity integer(pInt), intent(in) :: file - integer(pInt), parameter :: maxNchunks = 7_pInt - integer(pInt), dimension(1+2_pInt*maxNchunks) :: positions + integer(pInt), parameter :: maxNchunks = 7 + integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt) section, maxNinstance, i,j,k, output, mySize character(len=64) tag character(len=1024) line @@ -104,9 +103,9 @@ subroutine constitutive_j2_init(file) !$OMP END CRITICAL (write2out) maxNinstance = count(phase_constitution == constitutive_j2_label) - if (maxNinstance == 0_pInt) return + if (maxNinstance == 0) return - if (debug_verbosity > 0_pInt) then + if (debug_verbosity > 0) then !$OMP CRITICAL (write2out) write(6,'(a16,1x,i5)') '# instances:',maxNinstance write(6,*) @@ -127,12 +126,12 @@ subroutine constitutive_j2_init(file) allocate(constitutive_j2_n(maxNinstance)) ; constitutive_j2_n = 0.0_pReal allocate(constitutive_j2_h0(maxNinstance)) ; constitutive_j2_h0 = 0.0_pReal allocate(constitutive_j2_tausat(maxNinstance)) ; constitutive_j2_tausat = 0.0_pReal - allocate(constitutive_j2_a(maxNinstance)) ; constitutive_j2_a = 0.0_pReal + allocate(constitutive_j2_a(maxNinstance)) ; constitutive_j2_a = 0.0_pReal allocate(constitutive_j2_aTolResistance(maxNinstance)) ; constitutive_j2_aTolResistance = 0.0_pReal rewind(file) line = '' - section = 0_pInt + section = 0 do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to read(file,'(a1024)',END=100) line @@ -143,60 +142,62 @@ subroutine constitutive_j2_init(file) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1_pInt - output = 0_pInt ! reset output counter + section = section + 1 + output = 0 ! reset output counter endif - if (section > 0_pInt .and. phase_constitution(section) == constitutive_j2_label) then ! one of my sections + if (section > 0 .and. phase_constitution(section) == constitutive_j2_label) then ! one of my sections i = phase_constitutionInstance(section) ! which instance of my constitution is present phase positions = IO_stringPos(line,maxNchunks) tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key select case(tag) case ('(output)') - output = output + 1_pInt - constitutive_j2_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt)) + output = output + 1 + constitutive_j2_output(output,i) = IO_lc(IO_stringValue(line,positions,2)) case ('c11') - constitutive_j2_C11(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_C11(i) = IO_floatValue(line,positions,2) case ('c12') - constitutive_j2_C12(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_C12(i) = IO_floatValue(line,positions,2) case ('tau0') - constitutive_j2_tau0(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_tau0(i) = IO_floatValue(line,positions,2) case ('gdot0') - constitutive_j2_gdot0(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_gdot0(i) = IO_floatValue(line,positions,2) case ('n') - constitutive_j2_n(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_n(i) = IO_floatValue(line,positions,2) case ('h0') - constitutive_j2_h0(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_h0(i) = IO_floatValue(line,positions,2) case ('tausat') - constitutive_j2_tausat(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_tausat(i) = IO_floatValue(line,positions,2) case ('a', 'w0') - constitutive_j2_a(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_a(i) = IO_floatValue(line,positions,2) case ('taylorfactor') - constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2) case ('atol_resistance') - constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2_pInt) + constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2) + case default + call IO_error(210_pInt,ext_msg=tag) end select endif enddo -100 do i = 1_pInt,maxNinstance ! sanity checks - if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(210_pInt) - if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211_pInt) - if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(212_pInt) - if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(213_pInt) - if (constitutive_j2_a(i) <= 0.0_pReal) call IO_error(241_pInt) - if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(240_pInt) - if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(242_pInt) +100 do i = 1,maxNinstance ! sanity checks + if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(211_pInt,ext_msg='tau0') + if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='gdot0') + if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='n') + if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='tausat') + if (constitutive_j2_a(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='a') + if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='taylorfactor') + if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='aTol_resistance') enddo - do i = 1_pInt,maxNinstance - do j = 1_pInt,maxval(phase_Noutput) + do i = 1,maxNinstance + do j = 1,maxval(phase_Noutput) select case(constitutive_j2_output(j,i)) case('flowstress') - mySize = 1_pInt + mySize = 1_pInt case('strainrate') - mySize = 1_pInt + mySize = 1_pInt case default - mySize = 0_pInt + call IO_error(212_pInt,ext_msg=constitutive_j2_output(j,i)) end select if (mySize > 0_pInt) then ! any meaningful output found @@ -206,24 +207,22 @@ subroutine constitutive_j2_init(file) endif enddo - constitutive_j2_sizeDotState(i) = 1_pInt - constitutive_j2_sizeState(i) = 1_pInt + constitutive_j2_sizeDotState(i) = 1 + constitutive_j2_sizeState(i) = 1 - forall(k=1_pInt:3_pInt) - forall(j=1_pInt:3_pInt) & + forall(k=1:3) + forall(j=1:3) & constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i) constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i) - constitutive_j2_Cslip_66(k+3_pInt,k+3_pInt,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i)) + constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i)) end forall constitutive_j2_Cslip_66(1:6,1:6,i) = & math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(1:6,1:6,i))) - enddo return - endsubroutine @@ -244,7 +243,6 @@ pure function constitutive_j2_stateInit(myInstance) endfunction - !********************************************************************* !* relevant microstructural state * !********************************************************************* @@ -265,7 +263,6 @@ real(pReal), dimension(constitutive_j2_sizeState(myInstance)) :: & constitutive_j2_aTolState = constitutive_j2_aTolResistance(myInstance) - endfunction @@ -293,7 +290,6 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el) return - endfunction @@ -316,7 +312,6 @@ subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el) real(pReal) Temperature type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state - matID = phase_constitutionInstance(material_phase(ipc,ip,el)) endsubroutine @@ -342,7 +337,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v, material_phase, & phase_constitutionInstance - implicit none !*** input variables ***! @@ -466,7 +460,6 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el) endfunction - !**************************************************************** !* calculates the rate of change of temperature * !**************************************************************** diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index c6b87d0dd..6f136e867 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -341,13 +341,13 @@ do if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section section = section + 1_pInt - output = 0_pInt ! reset output counter + output = 0_pInt ! reset output counter cycle endif - if (section > 0_pInt .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections + if (section > 0_pInt .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections i = phase_constitutionInstance(section) ! which instance of my constitution is present phase positions = IO_stringPos(line,maxNchunks) - tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key select case(tag) case('constitution','/nonlocal/') cycle @@ -369,40 +369,40 @@ do case ('c44') constitutive_nonlocal_C44(i) = IO_floatValue(line,positions,2_pInt) case ('nslip') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_Nslip(f,i)& - = IO_intValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_Nslip(f,i) = IO_intValue(line,positions,1_pInt+f) case ('rhosgledgepos0') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgePos0(f,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_rhoSglEdgePos0(f,i) = IO_floatValue(line,positions,1_pInt+f) case ('rhosgledgeneg0') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgeNeg0(f,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_rhoSglEdgeNeg0(f,i) = IO_floatValue(line,positions,1_pInt+f) case ('rhosglscrewpos0') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewPos0(f,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_rhoSglScrewPos0(f,i) = IO_floatValue(line,positions,1_pInt+f) case ('rhosglscrewneg0') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewNeg0(f,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_rhoSglScrewNeg0(f,i) = IO_floatValue(line,positions,1_pInt+f) case ('rhodipedge0') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipEdge0(f,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_rhoDipEdge0(f,i) = IO_floatValue(line,positions,1_pInt+f) case ('rhodipscrew0') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipScrew0(f,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_rhoDipScrew0(f,i) = IO_floatValue(line,positions,1_pInt+f) case ('lambda0') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_lambda0PerSlipFamily(f,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_lambda0PerSlipFamily(f,i) = IO_floatValue(line,positions,1_pInt+f) case ('burgers') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_burgersPerSlipFamily(f,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_burgersPerSlipFamily(f,i) = IO_floatValue(line,positions,1_pInt+f) case('cutoffradius','r') constitutive_nonlocal_R(i) = IO_floatValue(line,positions,2_pInt) case('minimumdipoleheightedge','ddipminedge') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,1_pInt,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,1_pInt,i) = IO_floatValue(line,positions,1_pInt+f) case('minimumdipoleheightscrew','ddipminscrew') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,2_pInt,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,2_pInt,i) = IO_floatValue(line,positions,1_pInt+f) case('atomicvolume') constitutive_nonlocal_atomicVolume(i) = IO_floatValue(line,positions,2_pInt) case('selfdiffusionprefactor','dsd0') @@ -412,14 +412,14 @@ do case('atol_rho') constitutive_nonlocal_aTolRho(i) = IO_floatValue(line,positions,2_pInt) case ('interaction_slipslip') - forall (it = 1_pInt:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(it,i)& - = IO_floatValue(line,positions,1_pInt+it) + forall (it = 1_pInt:lattice_maxNinteraction) & + constitutive_nonlocal_interactionSlipSlip(it,i) = IO_floatValue(line,positions,1_pInt+it) case('peierlsstressedge') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_peierlsStressPerSlipFamily(f,1_pInt,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_peierlsStressPerSlipFamily(f,1_pInt,i) = IO_floatValue(line,positions,1_pInt+f) case('peierlsstressscrew') - forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_peierlsStressPerSlipFamily(f,2_pInt,i)& - = IO_floatValue(line,positions,1_pInt+f) + forall (f = 1_pInt:lattice_maxNslipFamily) & + constitutive_nonlocal_peierlsStressPerSlipFamily(f,2_pInt,i) = IO_floatValue(line,positions,1_pInt+f) case('doublekinkwidth') constitutive_nonlocal_doublekinkwidth(i) = IO_floatValue(line,positions,2_pInt) case('solidsolutionenergy') @@ -441,7 +441,7 @@ do case('surfacetransmissivity') constitutive_nonlocal_surfaceTransmissivity(i) = IO_floatValue(line,positions,2_pInt) case default - call IO_error(236_pInt,ext_msg=tag) + call IO_error(250_pInt,ext_msg=tag) end select endif enddo @@ -457,46 +457,46 @@ enddo !*** sanity checks if (myStructure < 1 .or. myStructure > 3) call IO_error(205_pInt) - if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(235_pInt,ext_msg='Nslip') + if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(251_pInt,ext_msg='Nslip') do o = 1,maxval(phase_Noutput) if(len(constitutive_nonlocal_output(o,i)) > 64) call IO_error(666_pInt) enddo do f = 1,lattice_maxNslipFamily if (constitutive_nonlocal_Nslip(f,i) > 0_pInt) then - if (constitutive_nonlocal_rhoSglEdgePos0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglEdgePos0') - if (constitutive_nonlocal_rhoSglEdgeNeg0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglEdgeNeg0') - if (constitutive_nonlocal_rhoSglScrewPos0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglScrewPos0') - if (constitutive_nonlocal_rhoSglScrewNeg0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglScrewNeg0') - if (constitutive_nonlocal_rhoDipEdge0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoDipEdge0') - if (constitutive_nonlocal_rhoDipScrew0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoDipScrew0') - if (constitutive_nonlocal_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='burgers') - if (constitutive_nonlocal_lambda0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='lambda0') + if (constitutive_nonlocal_rhoSglEdgePos0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglEdgePos0') + if (constitutive_nonlocal_rhoSglEdgeNeg0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglEdgeNeg0') + if (constitutive_nonlocal_rhoSglScrewPos0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglScrewPos0') + if (constitutive_nonlocal_rhoSglScrewNeg0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglScrewNeg0') + if (constitutive_nonlocal_rhoDipEdge0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoDipEdge0') + if (constitutive_nonlocal_rhoDipScrew0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoDipScrew0') + if (constitutive_nonlocal_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='burgers') + if (constitutive_nonlocal_lambda0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='lambda0') if (constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,1,i) <= 0.0_pReal) & - call IO_error(235_pInt,ext_msg='minimumDipoleHeightEdge') + call IO_error(251_pInt,ext_msg='minimumDipoleHeightEdge') if (constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,2,i) <= 0.0_pReal) & - call IO_error(235_pInt,ext_msg='minimumDipoleHeightScrew') - if (constitutive_nonlocal_peierlsStressPerSlipFamily(f,1,i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='peierlsStressEdge') - if (constitutive_nonlocal_peierlsStressPerSlipFamily(f,2,i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='peierlsStressScrew') + call IO_error(251_pInt,ext_msg='minimumDipoleHeightScrew') + if (constitutive_nonlocal_peierlsStressPerSlipFamily(f,1,i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='peierlsStressEdge') + if (constitutive_nonlocal_peierlsStressPerSlipFamily(f,2,i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='peierlsStressScrew') endif enddo if (any(constitutive_nonlocal_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 0.0_pReal)) & - call IO_error(235_pInt,ext_msg='interaction_SlipSlip') - if (constitutive_nonlocal_R(i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='r') - if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='atomicVolume') - if (constitutive_nonlocal_Dsd0(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='selfDiffusionPrefactor') - if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='selfDiffusionEnergy') - if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='aTol_rho') - if (constitutive_nonlocal_doublekinkwidth(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='doublekinkwidth') - if (constitutive_nonlocal_solidSolutionEnergy(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='solidSolutionEnergy') - if (constitutive_nonlocal_solidSolutionSize(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='solidSolutionSize') - if (constitutive_nonlocal_solidSolutionConcentration(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='solidSolutionConcentration') - if (constitutive_nonlocal_p(i) <= 0.0_pReal .or. constitutive_nonlocal_p(i) > 1.0_pReal) call IO_error(235_pInt,ext_msg='p') - if (constitutive_nonlocal_q(i) < 1.0_pReal .or. constitutive_nonlocal_q(i) > 2.0_pReal) call IO_error(235_pInt,ext_msg='q') - if (constitutive_nonlocal_viscosity(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='viscosity') - if (constitutive_nonlocal_fattack(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='attackFrequency') - if (constitutive_nonlocal_rhoSglScatter(i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglScatter') + call IO_error(251_pInt,ext_msg='interaction_SlipSlip') + if (constitutive_nonlocal_R(i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='r') + if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='atomicVolume') + if (constitutive_nonlocal_Dsd0(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='selfDiffusionPrefactor') + if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='selfDiffusionEnergy') + if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='aTol_rho') + if (constitutive_nonlocal_doublekinkwidth(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='doublekinkwidth') + if (constitutive_nonlocal_solidSolutionEnergy(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='solidSolutionEnergy') + if (constitutive_nonlocal_solidSolutionSize(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='solidSolutionSize') + if (constitutive_nonlocal_solidSolutionConcentration(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='solidSolutionConcentration') + if (constitutive_nonlocal_p(i) <= 0.0_pReal .or. constitutive_nonlocal_p(i) > 1.0_pReal) call IO_error(251_pInt,ext_msg='p') + if (constitutive_nonlocal_q(i) < 1.0_pReal .or. constitutive_nonlocal_q(i) > 2.0_pReal) call IO_error(251_pInt,ext_msg='q') + if (constitutive_nonlocal_viscosity(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='viscosity') + if (constitutive_nonlocal_fattack(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='attackFrequency') + if (constitutive_nonlocal_rhoSglScatter(i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglScatter') if (constitutive_nonlocal_surfaceTransmissivity(i) < 0.0_pReal & - .or. constitutive_nonlocal_surfaceTransmissivity(i) > 1.0_pReal) call IO_error(235_pInt,ext_msg='surfaceTransmissivity') + .or. constitutive_nonlocal_surfaceTransmissivity(i) > 1.0_pReal) call IO_error(251_pInt,ext_msg='surfaceTransmissivity') !*** determine total number of active slip systems @@ -649,7 +649,7 @@ do i = 1,maxNinstance case('dislocationstress') mySize = 6_pInt case default - call IO_error(237,ext_msg=constitutive_nonlocal_output(o,i)) + call IO_error(252_pInt,ext_msg=constitutive_nonlocal_output(o,i)) end select if (mySize > 0_pInt) then ! any meaningful output found @@ -662,13 +662,13 @@ do i = 1,maxNinstance !*** elasticity matrix and shear modulus according to material.config select case (myStructure) - case(1_pInt:2_pInt) ! cubic(s) + case(1_pInt:2_pInt) ! cubic(s) forall(k=1_pInt:3_pInt) forall(j=1_pInt:3_pInt) constitutive_nonlocal_Cslip_66(k,j,i) = constitutive_nonlocal_C12(i) constitutive_nonlocal_Cslip_66(k,k,i) = constitutive_nonlocal_C11(i) constitutive_nonlocal_Cslip_66(k+3_pInt,k+3_pInt,i) = constitutive_nonlocal_C44(i) end forall - case(3_pInt:) ! all hex + case(3_pInt:) ! all hex constitutive_nonlocal_Cslip_66(1,1,i) = constitutive_nonlocal_C11(i) constitutive_nonlocal_Cslip_66(2,2,i) = constitutive_nonlocal_C11(i) constitutive_nonlocal_Cslip_66(3,3,i) = constitutive_nonlocal_C33(i) @@ -1646,7 +1646,8 @@ logical considerEnteringFlux, & considerLeavingFlux #ifndef _OPENMP - if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then + if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) & + .or. .not. debug_selectiveDebugger)) then write(6,*) write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g write(6,*) @@ -1657,7 +1658,7 @@ select case(mesh_element(2,el)) case (1_pInt,6_pInt,7_pInt,8_pInt,9_pInt) ! all fine case default - call IO_error(-1_pInt,el,ip,g,'element type not supported for nonlocal constitution') + call IO_error(253_pInt,el,ip,g) end select myInstance = phase_constitutionInstance(material_phase(g,ip,el)) @@ -1699,11 +1700,12 @@ endif forall (t = 1_pInt:4_pInt) & gdot(1_pInt:ns,t) = rhoSgl(1_pInt:ns,t) * constitutive_nonlocal_burgers(1:ns,myInstance) * v(1:ns,t) -forall (s = 1_pInt:ns, t = 1_pInt:4_pInt, rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pReal) & ! contribution of used rho for changing sign of v +forall (s = 1_pInt:ns, t = 1_pInt:4_pInt, rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pReal) & ! contribution of used rho for changing sign of v gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgers(s,myInstance) * v(s,t) #ifndef _OPENMP - if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then + if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) & + .or. .not. debug_selectiveDebugger)) then write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot endif diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index d24960aea..f6d0b629d 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -345,6 +345,8 @@ subroutine constitutive_phenopowerlaw_init(file) case ('interaction_twintwin') forall (j = 1_pInt:lattice_maxNinteraction) & constitutive_phenopowerlaw_interaction_twintwin(j,i) = IO_floatValue(line,positions,1_pInt+j) + case default + call IO_error(220_pInt,ext_msg=tag) end select endif enddo @@ -362,21 +364,21 @@ subroutine constitutive_phenopowerlaw_init(file) constitutive_phenopowerlaw_totalNslip(i) = sum(constitutive_phenopowerlaw_Nslip(:,i)) ! how many slip systems altogether constitutive_phenopowerlaw_totalNtwin(i) = sum(constitutive_phenopowerlaw_Ntwin(:,i)) ! how many twin systems altogether - if (constitutive_phenopowerlaw_structure(i) < 1 ) call IO_error(205_pInt,i) + if (constitutive_phenopowerlaw_structure(i) < 1 ) call IO_error(205_pInt,e=i) if (any(constitutive_phenopowerlaw_tau0_slip(:,i) < 0.0_pReal .and. & - constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(210_pInt,i) - if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(211_pInt,i) - if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(212_pInt,i) + constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='tau0_slip') + if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(221_pInt,e=i,ext_msg='gdot0_slip') + if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(221_pInt,e=i,ext_msg='n_slip') if (any(constitutive_phenopowerlaw_tausat_slip(:,i) <= 0.0_pReal .and. & - constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(213_pInt,i) + constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='tausat_slip') if (any(constitutive_phenopowerlaw_a_slip(i) == 0.0_pReal .and. & - constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(214_pInt,i) + constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='a_slip') if (any(constitutive_phenopowerlaw_tau0_twin(:,i) < 0.0_pReal .and. & - constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(210_pInt,i) + constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='tau0_twin') if ( constitutive_phenopowerlaw_gdot0_twin(i) <= 0.0_pReal .and. & - any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211_pInt,i) + any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='gdot0_twin') if ( constitutive_phenopowerlaw_n_twin(i) <= 0.0_pReal .and. & - any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(212_pInt,i) + any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='n_twin') if (constitutive_phenopowerlaw_aTolResistance(i) <= 0.0_pReal) & constitutive_phenopowerlaw_aTolResistance(i) = 1.0_pReal ! default absolute tolerance 1 Pa @@ -418,7 +420,7 @@ subroutine constitutive_phenopowerlaw_init(file) ) mySize = 1_pInt case default - mySize = 0_pInt + call IO_error(222_pInt,ext_msg=constitutive_phenopowerlaw_output(j,i)) end select if (mySize > 0_pInt) then ! any meaningful output found diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index 3e2ca49e7..b3343b7ab 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -583,6 +583,8 @@ do ! read thru sections of forall (j = 1_pInt:lattice_maxNinteraction) & constitutive_titanmod_interactionTwinTwin(j,i) = IO_floatValue(line,positions,1_pInt+j) write(6,*) tag + case default + call IO_error(230_pInt,ext_msg=tag) end select endif enddo @@ -595,37 +597,38 @@ write(6,*) 'Material Property reading done' myStructure = constitutive_titanmod_structure(i) !* Sanity checks - if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt) - if (sum(constitutive_titanmod_Nslip(:,i)) <= 0_pInt) call IO_error(207_pInt) - if (sum(constitutive_titanmod_Ntwin(:,i)) < 0_pInt) call IO_error(208_pInt) !*** + if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt,e=i) + if (sum(constitutive_titanmod_Nslip(:,i)) <= 0_pInt) call IO_error(231_pInt,e=i,ext_msg='nslip') + if (sum(constitutive_titanmod_Ntwin(:,i)) < 0_pInt) call IO_error(231_pInt,e=i,ext_msg='ntwin') do f = 1,lattice_maxNslipFamily if (constitutive_titanmod_Nslip(f,i) > 0_pInt) then - if (constitutive_titanmod_rho_edge0(f,i) < 0.0_pReal) call IO_error(209_pInt) - if (constitutive_titanmod_rho_screw0(f,i) < 0.0_pReal) call IO_error(210_pInt) - if (constitutive_titanmod_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(211_pInt) - if (constitutive_titanmod_f0_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(212_pInt) - if (constitutive_titanmod_tau0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(229_pInt) - if (constitutive_titanmod_tau0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(233_pInt) - if (constitutive_titanmod_capre_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(234_pInt) - if (constitutive_titanmod_caprs_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235_pInt) - if (constitutive_titanmod_v0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt) - if (constitutive_titanmod_v0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt) - if (constitutive_titanmod_kinkcriticallength_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(238_pInt) + if (constitutive_titanmod_rho_edge0(f,i) < 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='rho_edge0') + if (constitutive_titanmod_rho_screw0(f,i) < 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='rho_screw0') + if (constitutive_titanmod_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='slipburgers') + if (constitutive_titanmod_f0_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='f0') + if (constitutive_titanmod_tau0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='tau0e') + if (constitutive_titanmod_tau0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='tau0s') + if (constitutive_titanmod_capre_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='capre') + if (constitutive_titanmod_caprs_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='carrs') + if (constitutive_titanmod_v0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='v0e') + if (constitutive_titanmod_v0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='v0s') + if (constitutive_titanmod_kinkcriticallength_PerSlipFamily(f,i) <= 0.0_pReal) & + call IO_error(231_pInt,e=i,ext_msg='kinkCriticalLength') endif enddo do f = 1,lattice_maxNtwinFamily if (constitutive_titanmod_Ntwin(f,i) > 0_pInt) then - if (constitutive_titanmod_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221_pInt) !*** - if (constitutive_titanmod_twinf0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(228_pInt) - if (constitutive_titanmod_twinshearconstant_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(228_pInt) - if (constitutive_titanmod_twintau0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(229_pInt) - if (constitutive_titanmod_twingamma0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt) + if (constitutive_titanmod_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinburgers') + if (constitutive_titanmod_twinf0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinf0') + if (constitutive_titanmod_twinshearconstant_PerTwinFamily(f,i) <= 0.0_pReal) & + call IO_error(231_pInt,e=i,ext_msg='twinshearconstant') + if (constitutive_titanmod_twintau0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twintau0') + if (constitutive_titanmod_twingamma0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twingamma0') endif enddo -! if (any(constitutive_titanmod_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) call IO_error(229) - if (constitutive_titanmod_dc(i) <= 0.0_pReal) call IO_error(231_pInt) - if (constitutive_titanmod_twinhpconstant(i) <= 0.0_pReal) call IO_error(232_pInt) - if (constitutive_titanmod_aTolRho(i) <= 0.0_pReal) call IO_error(233_pInt) + if (constitutive_titanmod_dc(i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='dc') + if (constitutive_titanmod_twinhpconstant(i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinhpconstant') + if (constitutive_titanmod_aTolRho(i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='aTolRho') !* Determine total number of active slip or twin systems constitutive_titanmod_Nslip(:,i) = min(lattice_NslipSystem(:,myStructure),constitutive_titanmod_Nslip(:,i)) @@ -637,7 +640,7 @@ enddo !* Allocation of variables whose size depends on the total number of active slip systems maxTotalNslip = maxval(constitutive_titanmod_totalNslip) -maxTotalNtwin = maxval(constitutive_titanmod_totalNtwin) +maxTotalNtwin = maxval(constitutive_titanmod_totalNtwin) write(6,*) 'maxTotalNslip',maxTotalNslip write(6,*) 'maxTotalNtwin',maxTotalNtwin allocate(constitutive_titanmod_burgersPerSlipSystem(maxTotalNslip, maxNinstance)) @@ -796,7 +799,7 @@ do i = 1_pInt,maxNinstance ) mySize = 1_pInt case default - mySize = 0_pInt + call IO_error(232_pInt,ext_msg=constitutive_titanmod_output(o,i)) end select if (mySize > 0_pInt) then ! any meaningful output found diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 5725aeeed..f5ee810e9 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -234,8 +234,8 @@ allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt -if (.not. IO_open_jobFile(file,material_localFileExt)) then ! no local material configuration present... - if (.not. IO_open_file(file,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file +if (.not. IO_open_jobFile_stat(file,material_localFileExt)) then ! no local material configuration present... + call IO_open_file(file,material_configFile) ! ...open material.config file endif line = '' section = 0_pInt @@ -299,7 +299,7 @@ enddo ! write description file for crystallite output -if(.not. IO_write_jobFile(file,'outputCrystallite')) call IO_error (50_pInt) ! problems in writing file +call IO_write_jobFile(file,'outputCrystallite') do p = 1_pInt,material_Ncrystallite write(file,*) diff --git a/code/debug.f90 b/code/debug.f90 index 7b319b2cb..3f3148b04 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -73,7 +73,7 @@ subroutine debug_init() nMPstate, & nHomog use IO, only: IO_error, & - IO_open_file, & + IO_open_file_stat, & IO_isBlank, & IO_stringPos, & IO_stringValue, & @@ -108,7 +108,7 @@ subroutine debug_init() allocate(debug_MaterialpointLoopDistribution(nHomog+1)) ; debug_MaterialpointLoopDistribution = 0_pInt ! try to open the config file - if(IO_open_file(fileunit,debug_configFile)) then + if(IO_open_file_stat(fileunit,debug_configFile)) then line = '' ! read variables from config file and overwrite parameters diff --git a/code/homogenization.f90 b/code/homogenization.f90 index dbd023641..d9455de1d 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -75,7 +75,7 @@ subroutine homogenization_init(Temperature) use prec, only: pReal,pInt use math, only: math_I3 use debug, only: debug_verbosity -use IO, only: IO_error, IO_open_file, IO_open_jobFile, IO_write_jobFile +use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips use material use constitutive, only: constitutive_maxSizePostResults @@ -95,8 +95,8 @@ logical knownHomogenization ! --- PARSE HOMOGENIZATIONS FROM CONFIG FILE --- -if (.not. IO_open_jobFile(fileunit,material_localFileExt)) then ! no local material configuration present... - if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100) ! ...and cannot open material.config file +if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present... + call IO_open_file(fileunit,material_configFile) ! ... open material.config file endif call homogenization_isostrain_init(fileunit) call homogenization_RGC_init(fileunit) @@ -105,9 +105,7 @@ close(fileunit) ! --- WRITE DESCRIPTION FILE FOR HOMOGENIZATION OUTPUT --- -if(.not. IO_write_jobFile(fileunit,'outputHomogenization')) then ! problems in writing file - call IO_error (50) -endif +call IO_write_jobFile(fileunit,'outputHomogenization') do p = 1,material_Nhomogenization i = homogenization_typeInstance(p) ! which instance of this homogenization type knownHomogenization = .true. ! assume valid @@ -189,7 +187,7 @@ end forall endif homogenization_sizePostResults(i,e) = homogenization_RGC_sizePostResults(myInstance) case default - call IO_error(201,ext_msg=homogenization_type(mesh_element(3,e))) ! unknown type 201 is homogenization! + call IO_error(500_pInt,ext_msg=homogenization_type(mesh_element(3,e))) ! unknown homogenization end select enddo enddo diff --git a/code/lattice.f90 b/code/lattice.f90 index 9945c9398..d4b1ecb32 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -728,7 +728,7 @@ subroutine lattice_init() !* Module initialization * !************************************** use, intrinsic :: iso_fortran_env - use IO, only: IO_open_file,IO_open_jobFile,IO_countSections,IO_countTagInPart,IO_error + use IO, only: IO_open_file,IO_open_jobFile_stat,IO_countSections,IO_countTagInPart,IO_error use material, only: material_configfile,material_localFileExt,material_partPhase use debug, only: debug_verbosity implicit none @@ -743,8 +743,8 @@ subroutine lattice_init() #include "compilation_info.f90" !$OMP END CRITICAL (write2out) - if (.not. IO_open_jobFile(fileunit,material_localFileExt)) then ! no local material configuration present... - if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file + if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present... + call IO_open_file(fileunit,material_configFile) ! ... open material.config file endif Nsections = IO_countSections(fileunit,material_partPhase) lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex @@ -922,7 +922,7 @@ function lattice_initializeStructure(struct,CoverA) if (processMe) then if (myStructure > lattice_Nstructure) & - call IO_error(666_pint,0_pInt,0_pInt,0_pInt,'structure index too large') ! check for memory leakage + call IO_error(666_pInt,0_pInt,0_pInt,0_pInt,'structure index too large') ! check for memory leakage do i = 1,myNslip ! store slip system vectors and Schmid matrix for my structure lattice_sd(1:3,i,myStructure) = sd(1:3,i) lattice_st(1:3,i,myStructure) = st(1:3,i) diff --git a/code/material.f90 b/code/material.f90 index a665a3354..30d5a874b 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -109,7 +109,7 @@ subroutine material_init() use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pReal,pInt - use IO, only: IO_error, IO_open_file, IO_open_jobFile + use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat use debug, only: debug_verbosity implicit none @@ -124,8 +124,8 @@ subroutine material_init() #include "compilation_info.f90" !$OMP END CRITICAL (write2out) - if (.not. IO_open_jobFile(fileunit,material_localFileExt)) then ! no local material configuration present... - if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file + if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present... + call IO_open_file(fileunit,material_configFile) ! ...open material.config file endif call material_parseHomogenization(fileunit,material_partHomogenization) if (debug_verbosity > 0) then @@ -163,16 +163,16 @@ subroutine material_init() if (microstructure_crystallite(i) < 1 .or. & microstructure_crystallite(i) > material_Ncrystallite) call IO_error(150_pInt,i) if (minval(microstructure_phase(1:microstructure_Nconstituents(i),i)) < 1 .or. & - maxval(microstructure_phase(1:microstructure_Nconstituents(i),i)) > material_Nphase) call IO_error(155_pInt,i) + maxval(microstructure_phase(1:microstructure_Nconstituents(i),i)) > material_Nphase) call IO_error(151_pInt,i) if (minval(microstructure_texture(1:microstructure_Nconstituents(i),i)) < 1 .or. & - maxval(microstructure_texture(1:microstructure_Nconstituents(i),i)) > material_Ntexture) call IO_error(160_pInt,i) + maxval(microstructure_texture(1:microstructure_Nconstituents(i),i)) > material_Ntexture) call IO_error(152_pInt,i) if (abs(sum(microstructure_fraction(:,i)) - 1.0_pReal) >= 1.0e-10_pReal) then if (debug_verbosity > 0) then !$OMP CRITICAL (write2out) write(6,*)'sum of microstructure fraction = ',sum(microstructure_fraction(:,i)) !$OMP END CRITICAL (write2out) endif - call IO_error(170_pInt,i) + call IO_error(153_pInt,i) endif enddo if (debug_verbosity > 0) then @@ -227,7 +227,7 @@ subroutine material_parseHomogenization(file,myPart) Nsections = IO_countSections(file,myPart) material_Nhomogenization = Nsections - if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart) + if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) allocate(homogenization_name(Nsections)); homogenization_name = '' allocate(homogenization_type(Nsections)); homogenization_type = '' @@ -295,7 +295,7 @@ subroutine material_parseMicrostructure(file,myPart) Nsections = IO_countSections(file,myPart) material_Nmicrostructure = Nsections - if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart) + if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) allocate(microstructure_name(Nsections)); microstructure_name = '' allocate(microstructure_crystallite(Nsections)); microstructure_crystallite = 0_pInt @@ -372,7 +372,7 @@ subroutine material_parseCrystallite(file,myPart) Nsections = IO_countSections(file,myPart) material_Ncrystallite = Nsections - if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart) + if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) allocate(crystallite_name(Nsections)); crystallite_name = '' allocate(crystallite_Noutput(Nsections)); crystallite_Noutput = 0_pInt @@ -418,7 +418,7 @@ subroutine material_parsePhase(file,myPart) Nsections = IO_countSections(file,myPart) material_Nphase = Nsections - if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart) + if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) allocate(phase_name(Nsections)); phase_name = '' allocate(phase_constitution(Nsections)); phase_constitution = '' @@ -482,7 +482,7 @@ subroutine material_parseTexture(file,myPart) Nsections = IO_countSections(file,myPart) material_Ntexture = Nsections - if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart) + if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) allocate(texture_name(Nsections)); texture_name = '' allocate(texture_ODFfile(Nsections)); texture_ODFfile = '' @@ -645,9 +645,9 @@ subroutine material_populateGrains() homog = mesh_element(3,e) micro = mesh_element(4,e) if (homog < 1 .or. homog > material_Nhomogenization) & ! out of bounds - call IO_error(130_pInt,e,0_pInt,0_pInt) + call IO_error(154_pInt,e,0_pInt,0_pInt) if (micro < 1 .or. micro > material_Nmicrostructure) & ! out of bounds - call IO_error(140_pInt,e,0_pInt,0_pInt) + call IO_error(155_pInt,e,0_pInt,0_pInt) if (microstructure_elemhomo(micro)) then dGrains = homogenization_Ngrains(homog) else @@ -762,7 +762,7 @@ subroutine material_populateGrains() else ! hybrid IA ! --------- orientationOfGrain(:,grain+1:grain+myNorientations) = IO_hybridIA(myNorientations,texture_ODFfile(textureID)) - if (all(orientationOfGrain(:,grain+1) == -1.0_pReal)) call IO_error(105_pInt) + if (all(orientationOfGrain(:,grain+1) == -1.0_pReal)) call IO_error(156_pInt) constituentGrain = constituentGrain + myNorientations endif diff --git a/code/math.f90 b/code/math.f90 index 75b2212f2..5201c04bf 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -201,7 +201,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = & if ( any(abs( q-q2) > tol_math_check) .and. & any(abs(-q-q2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(670_pInt,ext_msg=error_msg) + call IO_error(401_pInt,ext_msg=error_msg) endif ! +++ q -> R -> q +++ @@ -210,7 +210,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = & if ( any(abs( q-q2) > tol_math_check) .and. & any(abs(-q-q2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(671_pInt,ext_msg=error_msg) + call IO_error(402_pInt,ext_msg=error_msg) endif ! +++ q -> euler -> q +++ @@ -219,7 +219,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = & if ( any(abs( q-q2) > tol_math_check) .and. & any(abs(-q-q2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(672_pInt,ext_msg=error_msg) + call IO_error(403_pInt,ext_msg=error_msg) endif ! +++ R -> euler -> R +++ @@ -227,7 +227,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = & R2 = math_EulerToR(Eulers) if ( any(abs( R-R2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) 'maximum deviation ',maxval(abs( R-R2)) - call IO_error(673_pInt,ext_msg=error_msg) + call IO_error(404_pInt,ext_msg=error_msg) endif ENDSUBROUTINE math_init @@ -1820,7 +1820,7 @@ function math_QuaternionDisorientation(Q1, Q2, symmetryType) enddo; enddo; enddo case default - call IO_error(550_pInt,symmetryType) ! complain about unknown symmetry + call IO_error(450_pInt,symmetryType) ! complain about unknown symmetry end select endfunction math_QuaternionDisorientation @@ -2553,7 +2553,7 @@ end subroutine r(1:ndim) = 0.0_pReal - if (any (base(1:ndim) <= 1_pInt)) call IO_error(error_ID=801_pInt) + if (any (base(1:ndim) <= 1_pInt)) call IO_error(error_ID=405_pInt) base_inv(1:ndim) = 1.0_pReal / real (base(1:ndim), pReal) @@ -2795,7 +2795,7 @@ end subroutine else if (n <= prime_max) then prime = npvec(n) else - call IO_error(error_ID=802_pInt) + call IO_error(error_ID=406_pInt) end if endfunction prime @@ -3262,7 +3262,7 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) step = geomdim/real(res, pReal) - if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt) + if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) call fftw_set_timelimit(fftw_timelimit) defgrad_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*9_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) call c_f_pointer(defgrad_fftw, defgrad_real, [res(1)+2_pInt,res(2),res(3),3_pInt,3_pInt]) @@ -3385,7 +3385,7 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl) wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) - if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt) + if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) call fftw_set_timelimit(fftw_timelimit) field_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) @@ -3500,10 +3500,10 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence) print '(a,3(i5))', ' Resolution:', res endif - res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) + res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal) -if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt) +if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) call fftw_set_timelimit(fftw_timelimit) field_fftw = fftw_alloc_complex(int(res1_red*res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) @@ -3512,9 +3512,9 @@ if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt) call c_f_pointer(divergence_fftw, divergence_real, [res(1)+2_pInt,res(2),res(3),vec_tens]) call c_f_pointer(divergence_fftw, divergence_fourier,[res1_red ,res(2),res(3),vec_tens]) - fftw_forth = fftw_plan_many_dft_r2c(3_pInt,(/res(3),res(2) ,res(1)/),vec_tens*3_pInt,& ! dimensions , length in each dimension in reversed order + fftw_forth = fftw_plan_many_dft_r2c(3_pInt,(/res(3),res(2) ,res(1)/),vec_tens*3_pInt,& ! dimensions , length in each dimension in reversed order field_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order - 1_pInt, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions + 1_pInt, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions field_fourier,(/res(3),res(2) ,res1_red/),& 1_pInt, res(3)*res(2)* res1_red,fftw_planner_flag) @@ -3522,7 +3522,7 @@ if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt) divergence_fourier,(/res(3),res(2) ,res1_red/),& 1_pInt, res(3)*res(2)* res1_red,& divergence_real,(/res(3),res(2) ,res(1)+2_pInt/),& - 1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) ! padding + 1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) ! padding do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) field_real(i,j,k,1:vec_tens,1:3) = field(i,j,k,1:vec_tens,1:3) ! ensure that data is aligned properly (fftw_alloc) enddo; enddo; enddo @@ -3559,12 +3559,12 @@ if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt) call fftw_execute_dft_c2r(fftw_back, divergence_fourier, divergence_real) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - divergence(i,j,k,1:vec_tens) = divergence_real(i,j,k,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc) + divergence(i,j,k,1:vec_tens) = divergence_real(i,j,k,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc) enddo; enddo; enddo divergence = divergence * wgt call fftw_destroy_plan(fftw_forth); call fftw_destroy_plan(fftw_back) - call c_f_pointer(C_NULL_PTR, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) ! let all pointers point on NULL-Type + call c_f_pointer(C_NULL_PTR, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) ! let all pointers point on NULL-Type call c_f_pointer(C_NULL_PTR, field_fourier, [res1_red ,res(2),res(3),vec_tens,3_pInt]) call c_f_pointer(C_NULL_PTR, divergence_real, [res(1)+2_pInt,res(2),res(3),vec_tens]) call c_f_pointer(C_NULL_PTR, divergence_fourier,[res1_red ,res(2),res(3),vec_tens]) diff --git a/code/mesh.f90 b/code/mesh.f90 index 427e29438..a69a1e858 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -281,59 +281,56 @@ call mesh_build_FEdata() ! --- get properties of the different types of elements - if (IO_open_inputFile(fileUnit,FEmodelGeometry)) then ! --- parse info from input file... + call IO_open_inputFile(fileUnit,FEmodelGeometry) ! --- parse info from input file... - select case (FEsolver) - case ('Spectral') - call mesh_spectral_count_nodesAndElements(fileUnit) - call mesh_spectral_count_cpElements() - call mesh_spectral_map_elements() - call mesh_spectral_map_nodes() - call mesh_spectral_count_cpSizes() - call mesh_spectral_build_nodes(fileUnit) - call mesh_spectral_build_elements(fileUnit) - - case ('Marc') - call mesh_marc_get_tableStyles(fileUnit) - call mesh_marc_count_nodesAndElements(fileUnit) - call mesh_marc_count_elementSets(fileUnit) - call mesh_marc_map_elementSets(fileUnit) - call mesh_marc_count_cpElements(fileUnit) - call mesh_marc_map_elements(fileUnit) - call mesh_marc_map_nodes(fileUnit) - call mesh_marc_build_nodes(fileUnit) - call mesh_marc_count_cpSizes(fileunit) - call mesh_marc_build_elements(fileUnit) - case ('Abaqus') - noPart = IO_abaqus_hasNoPart(fileUnit) - call mesh_abaqus_count_nodesAndElements(fileUnit) - call mesh_abaqus_count_elementSets(fileUnit) - call mesh_abaqus_count_materials(fileUnit) - call mesh_abaqus_map_elementSets(fileUnit) - call mesh_abaqus_map_materials(fileUnit) - call mesh_abaqus_count_cpElements(fileUnit) - call mesh_abaqus_map_elements(fileUnit) - call mesh_abaqus_map_nodes(fileUnit) - call mesh_abaqus_build_nodes(fileUnit) - call mesh_abaqus_count_cpSizes(fileunit) - call mesh_abaqus_build_elements(fileUnit) - end select - call mesh_get_damaskOptions(fileUnit) - close (fileUnit) + select case (FEsolver) + case ('Spectral') + call mesh_spectral_count_nodesAndElements(fileUnit) + call mesh_spectral_count_cpElements() + call mesh_spectral_map_elements() + call mesh_spectral_map_nodes() + call mesh_spectral_count_cpSizes() + call mesh_spectral_build_nodes(fileUnit) + call mesh_spectral_build_elements(fileUnit) - call mesh_build_subNodeCoords() - call mesh_build_ipCoordinates() - call mesh_build_ipVolumes() - call mesh_build_ipAreas() - call mesh_build_nodeTwins() - call mesh_build_sharedElems() - call mesh_build_ipNeighborhood() - call mesh_tell_statistics() + case ('Marc') + call mesh_marc_get_tableStyles(fileUnit) + call mesh_marc_count_nodesAndElements(fileUnit) + call mesh_marc_count_elementSets(fileUnit) + call mesh_marc_map_elementSets(fileUnit) + call mesh_marc_count_cpElements(fileUnit) + call mesh_marc_map_elements(fileUnit) + call mesh_marc_map_nodes(fileUnit) + call mesh_marc_build_nodes(fileUnit) + call mesh_marc_count_cpSizes(fileunit) + call mesh_marc_build_elements(fileUnit) + case ('Abaqus') + noPart = IO_abaqus_hasNoPart(fileUnit) + call mesh_abaqus_count_nodesAndElements(fileUnit) + call mesh_abaqus_count_elementSets(fileUnit) + call mesh_abaqus_count_materials(fileUnit) + call mesh_abaqus_map_elementSets(fileUnit) + call mesh_abaqus_map_materials(fileUnit) + call mesh_abaqus_count_cpElements(fileUnit) + call mesh_abaqus_map_elements(fileUnit) + call mesh_abaqus_map_nodes(fileUnit) + call mesh_abaqus_build_nodes(fileUnit) + call mesh_abaqus_count_cpSizes(fileunit) + call mesh_abaqus_build_elements(fileUnit) + end select + call mesh_get_damaskOptions(fileUnit) + close (fileUnit) + + call mesh_build_subNodeCoords() + call mesh_build_ipCoordinates() + call mesh_build_ipVolumes() + call mesh_build_ipAreas() + call mesh_build_nodeTwins() + call mesh_build_sharedElems() + call mesh_build_ipNeighborhood() + call mesh_tell_statistics() - parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive - else - call IO_error(error_ID=101_pInt) ! cannot open input file - endif + parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive FEsolving_execElem = [ 1_pInt,mesh_NcpElems] allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt @@ -1487,7 +1484,7 @@ enddo if (keyword(1:4) == 'head') then headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt else - call IO_error(error_ID=42_pInt) + call IO_error(error_ID=842_pInt) endif rewind(myUnit) @@ -2428,7 +2425,7 @@ subroutine mesh_marc_count_cpSizes (myUnit) if (keyword(1:4) == 'head') then headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt else - call IO_error(error_ID=42_pInt) + call IO_error(error_ID=842_pInt) endif rewind(myUnit) @@ -2465,9 +2462,9 @@ subroutine mesh_marc_count_cpSizes (myUnit) ! --- sanity checks --- - if ((.not. gotDimension) .or. (.not. gotResolution)) call IO_error(error_ID=42_pInt) - if ((a < 1) .or. (b < 1) .or. (c < 0)) call IO_error(error_ID=43_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=44_pInt) + if ((.not. gotDimension) .or. (.not. gotResolution)) call IO_error(error_ID=842_pInt) + if ((a < 1) .or. (b < 1) .or. (c < 0)) 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) forall (n = 0:mesh_Nnodes-1) mesh_node0(1,n+1) = x * dble(mod(n,a)) / dble(a-1_pInt) @@ -2616,7 +2613,7 @@ subroutine mesh_marc_count_cpSizes (myUnit) if (keyword(1:4) == 'head') then headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt else - call IO_error(error_ID=42_pInt) + call IO_error(error_ID=842_pInt) endif rewind(myUnit) @@ -2679,7 +2676,7 @@ subroutine mesh_marc_count_cpSizes (myUnit) enddo 110 deallocate(microstructures) - if (e /= mesh_NcpElems) call IO_error(180_pInt,e) + if (e /= mesh_NcpElems) call IO_error(880_pInt,e) endsubroutine @@ -3365,13 +3362,13 @@ character(len=64) fmt integer(pInt) i,e,n,f,t -if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=110_pInt) ! no homogenization specified -if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=120_pInt) ! no microstructure specified +if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=170_pInt) ! no homogenization specified +if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=180_pInt) ! no microstructure specified allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt do e = 1,mesh_NcpElems - if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=110_pInt,e=e) ! no homogenization specified - if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=120_pInt,e=e) ! no microstructure specified + if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=170_pInt,e=e) ! no homogenization specified + if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=180_pInt,e=e) ! no microstructure specified mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = & mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1 ! count combinations of homogenization and microstructure enddo @@ -3490,4 +3487,4 @@ deallocate(mesh_HomogMicro) endsubroutine -END MODULE mesh \ No newline at end of file +END MODULE mesh diff --git a/code/numerics.f90 b/code/numerics.f90 index b38f2e018..786ff4c82 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -97,7 +97,7 @@ subroutine numerics_init() use prec, only: pInt, & pReal use IO, only: IO_error, & - IO_open_file, & + IO_open_file_stat, & IO_isBlank, & IO_stringPos, & IO_stringValue, & @@ -133,15 +133,17 @@ subroutine numerics_init() !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! ...and use it as number of threads for parallel execution ! try to open the config file - if(IO_open_file(fileunit,numerics_configFile)) then + if(IO_open_file_stat(fileunit,numerics_configFile)) then !$OMP CRITICAL (write2out) write(6,*) ' ... using values from config file' write(6,*) !$OMP END CRITICAL (write2out) + + !* read variables from config file and overwrite parameters + line = '' - ! read variables from config file and overwrite parameters do read(fileunit,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines @@ -195,7 +197,8 @@ subroutine numerics_init() case ('integratorstiffness') numerics_integrator(2) = IO_intValue(line,positions,2_pInt) -!* RGC parameters: + !* RGC parameters: + case ('atol_rgc') absTol_RGC = IO_floatValue(line,positions,2_pInt) case ('rtol_rgc') @@ -223,7 +226,8 @@ subroutine numerics_init() case ('discrepancypower_rgc') volDiscrPow_RGC = IO_floatValue(line,positions,2_pInt) -!* spectral parameters + !* spectral parameters + case ('err_div_tol') err_div_tol = IO_floatValue(line,positions,2_pInt) case ('err_stress_tolrel') @@ -247,9 +251,13 @@ subroutine numerics_init() case ('cut_off_value') cut_off_value = IO_floatValue(line,positions,2_pInt) -!* Random seeding parameters + !* Random seeding parameters + case ('fixed_seed') fixedSeed = IO_intValue(line,positions,2_pInt) + + case default + call IO_error(300_pInt,ext_msg=tag) endselect enddo 100 close(fileunit) @@ -263,6 +271,7 @@ subroutine numerics_init() !$OMP END CRITICAL (write2out) endif + select case(IO_lc(fftw_planner_string)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution fftw_planner_flag = 64_pInt @@ -277,8 +286,11 @@ subroutine numerics_init() fftw_planner_flag = 32_pInt end select - ! writing parameters to output file + + !* writing parameters to output file + !$OMP CRITICAL (write2out) + write(6,'(a24,1x,e8.1)') ' relevantStrain: ',relevantStrain write(6,'(a24,1x,e8.1)') ' defgradTolerance: ',defgradTolerance write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness @@ -303,7 +315,8 @@ subroutine numerics_init() write(6,'(a24,1x,e8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog write(6,'(a24,1x,i8,/)') ' nMPstate: ',nMPstate -!* RGC parameters + !* RGC parameters + write(6,'(a24,1x,e8.1)') ' aTol_RGC: ',absTol_RGC write(6,'(a24,1x,e8.1)') ' rTol_RGC: ',relTol_RGC write(6,'(a24,1x,e8.1)') ' aMax_RGC: ',absMax_RGC @@ -317,7 +330,8 @@ subroutine numerics_init() write(6,'(a24,1x,e8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC write(6,'(a24,1x,e8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC -!* spectral parameters + !* spectral parameters + write(6,'(a24,1x,e8.1)') ' err_div_tol: ',err_div_tol write(6,'(a24,1x,e8.1)') ' err_stress_tolrel: ',err_stress_tolrel write(6,'(a24,1x,i8)') ' itmax: ',itmax @@ -334,58 +348,66 @@ subroutine numerics_init() write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma write(6,'(a24,1x,L8,/)') ' simplified_algorithm: ',simplified_algorithm write(6,'(a24,1x,e8.1)') ' cut_off_value: ',cut_off_value -!* Random seeding parameters + + !* Random seeding parameters + write(6,'(a24,1x,i16,/)') ' fixed_seed: ',fixedSeed + !$OMP END CRITICAL (write2out) -!* openMP parameter + + !* openMP parameter !$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt - ! sanity check - if (relevantStrain <= 0.0_pReal) call IO_error(260_pInt) - if (defgradTolerance <= 0.0_pReal) call IO_error(294_pInt) - if (iJacoStiffness < 1_pInt) call IO_error(261_pInt) - if (iJacoLpresiduum < 1_pInt) call IO_error(262_pInt) - if (pert_Fg <= 0.0_pReal) call IO_error(263_pInt) + + !* sanity check + + if (relevantStrain <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relevantStrain') + if (defgradTolerance <= 0.0_pReal) call IO_error(301_pInt,ext_msg='defgradTolerance') + if (iJacoStiffness < 1_pInt) call IO_error(301_pInt,ext_msg='iJacoStiffness') + if (iJacoLpresiduum < 1_pInt) call IO_error(301_pInt,ext_msg='iJacoLpresiduum') + if (pert_Fg <= 0.0_pReal) call IO_error(301_pInt,ext_msg='pert_Fg') if (pert_method <= 0_pInt .or. pert_method >= 4_pInt) & - call IO_error(299_pInt) - if (nHomog < 1_pInt) call IO_error(264_pInt) - if (nMPstate < 1_pInt) call IO_error(279_pInt) !! missing in IO !! - if (nCryst < 1_pInt) call IO_error(265_pInt) - if (nState < 1_pInt) call IO_error(266_pInt) - if (nStress < 1_pInt) call IO_error(267_pInt) - if (subStepMinCryst <= 0.0_pReal) call IO_error(268_pInt) - if (subStepSizeCryst <= 0.0_pReal) call IO_error(268_pInt) - if (stepIncreaseCryst <= 0.0_pReal) call IO_error(268_pInt) - if (subStepMinHomog <= 0.0_pReal) call IO_error(268_pInt) - if (subStepSizeHomog <= 0.0_pReal) call IO_error(268_pInt) - if (stepIncreaseHomog <= 0.0_pReal) call IO_error(268_pInt) - if (rTol_crystalliteState <= 0.0_pReal) call IO_error(269_pInt) - if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276_pInt) !! oops !! - if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270_pInt) - if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271_pInt) + call IO_error(301_pInt,ext_msg='pert_method') + if (nHomog < 1_pInt) call IO_error(301_pInt,ext_msg='nHomog') + if (nMPstate < 1_pInt) call IO_error(301_pInt,ext_msg='nMPstate') + if (nCryst < 1_pInt) call IO_error(301_pInt,ext_msg='nCryst') + if (nState < 1_pInt) call IO_error(301_pInt,ext_msg='nState') + if (nStress < 1_pInt) call IO_error(301_pInt,ext_msg='nStress') + if (subStepMinCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinCryst') + if (subStepSizeCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeCryst') + if (stepIncreaseCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseCryst') + if (subStepMinHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinHomog') + if (subStepSizeHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeHomog') + if (stepIncreaseHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseHomog') + if (rTol_crystalliteState <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteState') + if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteTemperature') + if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteStress') + if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='aTol_crystalliteStress') if (any(numerics_integrator <= 0_pInt) .or. any(numerics_integrator >= 6_pInt)) & - call IO_error(298_pInt) + call IO_error(301_pInt,ext_msg='integrator') + !* RGC parameters - if (absTol_RGC <= 0.0_pReal) call IO_error(272_pInt) - if (relTol_RGC <= 0.0_pReal) call IO_error(273_pInt) - if (absMax_RGC <= 0.0_pReal) call IO_error(274_pInt) - if (relMax_RGC <= 0.0_pReal) call IO_error(275_pInt) - if (pPert_RGC <= 0.0_pReal) call IO_error(276_pInt) !! oops !! - if (xSmoo_RGC <= 0.0_pReal) call IO_error(277_pInt) - if (viscPower_RGC < 0.0_pReal) call IO_error(278_pInt) - if (viscModus_RGC < 0.0_pReal) call IO_error(278_pInt) - if (refRelaxRate_RGC <= 0.0_pReal) call IO_error(278_pInt) - if (maxdRelax_RGC <= 0.0_pReal) call IO_error(288_pInt) - if (maxVolDiscr_RGC <= 0.0_pReal) call IO_error(289_pInt) - if (volDiscrMod_RGC < 0.0_pReal) call IO_error(289_pInt) - if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(289_pInt) + if (absTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absTol_RGC') + if (relTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relTol_RGC') + if (absMax_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absMax_RGC') + if (relMax_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relMax_RGC') + if (pPert_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='pPert_RGC') + if (xSmoo_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='xSmoo_RGC') + if (viscPower_RGC < 0.0_pReal) call IO_error(301_pInt,ext_msg='viscPower_RGC') + if (viscModus_RGC < 0.0_pReal) call IO_error(301_pInt,ext_msg='viscModus_RGC') + if (refRelaxRate_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='refRelaxRate_RGC') + if (maxdRelax_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='maxdRelax_RGC') + if (maxVolDiscr_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='maxVolDiscr_RGC') + if (volDiscrMod_RGC < 0.0_pReal) call IO_error(301_pInt,ext_msg='volDiscrMod_RGC') + if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='volDiscrPw_RGC') -!* spectral parameters - if (err_div_tol <= 0.0_pReal) call IO_error(49_pInt) - if (err_stress_tolrel <= 0.0_pReal) call IO_error(49_pInt) - if (itmax <= 1.0_pInt) call IO_error(49_pInt) + !* spectral parameters + + if (err_div_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_div_tol') + if (err_stress_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolrel') + if (itmax <= 1.0_pInt) call IO_error(301_pInt,ext_msg='itmax') if (fixedSeed <= 0_pInt) then !$OMP CRITICAL (write2out)