From 3cb279b083efe798bacffb0e648d7bedc3996b21 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 31 Aug 2018 09:39:33 +0200 Subject: [PATCH] one more loop not needed --- src/DAMASK_spectral.f90 | 158 +++++++++++++++++++--------------------- 1 file changed, 74 insertions(+), 84 deletions(-) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index 92012dcfa..fb67d575d 100644 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -115,7 +115,7 @@ program DAMASK_spectral stagIterate integer(pInt) :: & i, j, k, l, field, & - errorID, & + errorID = 0_pInt, & cutBackLevel = 0_pInt, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ stepFraction = 0_pInt !< fraction of current time interval integer(pInt) :: & @@ -166,6 +166,7 @@ program DAMASK_spectral if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1 if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 allocate(solres(nActiveFields)) + allocate(newLoadCase%ID(nActiveFields)) !-------------------------------------------------------------------------------------------------- ! assign mechanics solver depending on selected type @@ -188,7 +189,7 @@ program DAMASK_spectral end select !-------------------------------------------------------------------------------------------------- -! reading basic information from load case file and allocate data structure containing load cases +! reading information from load case file and to sanity checks allocate (loadCases(0)) ! array of load cases open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read') if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=trim(loadCaseFile)) @@ -209,19 +210,17 @@ program DAMASK_spectral enddo ! count all identifiers to allocate memory and do sanity check if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1_pInt) & ! sanity check call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase - loadCases = [loadCases,newLoadCase] - currentLoadCase = currentLoadCase + 1_pInt - loadCases(currentLoadCase)%stress%myType='stress' - allocate(loadCases(currentLoadCase)%ID(nActiveFields)) + + newLoadCase%stress%myType='stress' field = 1 - loadCases(currentLoadCase)%ID(field) = FIELD_MECH_ID ! mechanical active by default + newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then field = field + 1 - loadCases(currentLoadCase)%ID(field) = FIELD_THERMAL_ID + newLoadCase%ID(field) = FIELD_THERMAL_ID endif thermalActive damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then field = field + 1 - loadCases(currentLoadCase)%ID(field) = FIELD_DAMAGE_ID + newLoadCase%ID(field) = FIELD_DAMAGE_ID endif damageActive do i = 1_pInt, chunkPos(1) @@ -230,46 +229,43 @@ program DAMASK_spectral temp_valueVector = 0.0_pReal if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'fdot'.or. & ! in case of Fdot, set type to fdot IO_lc(IO_stringValue(line,chunkPos,i)) == 'dotf') then - loadCases(currentLoadCase)%deformation%myType = 'fdot' + newLoadCase%deformation%myType = 'fdot' else if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'f') then - loadCases(currentLoadCase)%deformation%myType = 'f' + newLoadCase%deformation%myType = 'f' else - loadCases(currentLoadCase)%deformation%myType = 'l' + newLoadCase%deformation%myType = 'l' endif do j = 1_pInt, 9_pInt temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a * if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo - loadCases(currentLoadCase)%deformation%maskLogical = & ! logical mask in 3x3 notation - transpose(reshape(temp_maskVector,[ 3,3])) - loadCases(currentLoadCase)%deformation%maskFloat = & ! float (1.0/0.0) mask in 3x3 notation - merge(ones,zeros,loadCases(currentLoadCase)%deformation%maskLogical) - loadCases(currentLoadCase)%deformation%values = math_plain9to33(temp_valueVector) ! values in 3x3 notation + newLoadCase%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) ! logical mask in 3x3 notation + newLoadCase%deformation%maskFloat = merge(ones,zeros,newLoadCase%deformation%maskLogical)! float (1.0/0.0) mask in 3x3 notation + newLoadCase%deformation%values = math_plain9to33(temp_valueVector) ! values in 3x3 notation case('p','pk1','piolakirchhoff','stress', 's') temp_valueVector = 0.0_pReal do j = 1_pInt, 9_pInt temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo - loadCases(currentLoadCase)%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) - loadCases(currentLoadCase)%stress%maskFloat = merge(ones,zeros,& - loadCases(currentLoadCase)%stress%maskLogical) - loadCases(currentLoadCase)%stress%values = math_plain9to33(temp_valueVector) + newLoadCase%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) + newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical) + newLoadCase%stress%values = math_plain9to33(temp_valueVector) case('t','time','delta') ! increment time - loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1_pInt) + newLoadCase%time = IO_floatValue(line,chunkPos,i+1_pInt) case('n','incs','increments','steps') ! number of increments - loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt) + newLoadCase%incs = IO_intValue(line,chunkPos,i+1_pInt) case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) - loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt) - loadCases(currentLoadCase)%logscale = 1_pInt + newLoadCase%incs = IO_intValue(line,chunkPos,i+1_pInt) + newLoadCase%logscale = 1_pInt case('freq','frequency','outputfreq') ! frequency of result writings - loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt) + newLoadCase%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt) case('r','restart','restartwrite') ! frequency of writing restart information - loadCases(currentLoadCase)%restartfrequency = & + newLoadCase%restartfrequency = & max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt)) case('guessreset','dropguessing') - loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory - case('euler') ! rotation of currentLoadCase given in euler angles + newLoadCase%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory + case('euler') ! rotation of load case given in euler angles temp_valueVector = 0.0_pReal l = 1_pInt ! assuming values given in degrees k = 1_pInt ! assuming keyword indicating degree/radians present @@ -284,87 +280,79 @@ program DAMASK_spectral temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j) enddo if (l == 1_pInt) temp_valueVector(1:3) = temp_valueVector(1:3) * inRad ! convert to rad - loadCases(currentLoadCase)%rotation = math_EulerToR(temp_valueVector(1:3)) ! convert rad Eulers to rotation matrix - case('rotation','rot') ! assign values for the rotation of currentLoadCase matrix + newLoadCase%rotation = math_EulerToR(temp_valueVector(1:3)) ! convert rad Eulers to rotation matrix + case('rotation','rot') ! assign values for the rotation matrix temp_valueVector = 0.0_pReal do j = 1_pInt, 9_pInt temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) enddo - loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) + newLoadCase%rotation = math_plain9to33(temp_valueVector) end select enddo - enddo + currentLoadCase = currentLoadCase + 1_pInt + if(currentLoadCase == 1_pInt) newLoadCase%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first load case - close(fileUnit) - -!-------------------------------------------------------------------------------------------------- -! consistency checks and output of load case - loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase - errorID = 0_pInt - if (worldrank == 0) then - checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases) + reportAndCheck: if (worldrank == 0) then write (loadcase_string, '(i6)' ) currentLoadCase write(6,'(1x,a,i6)') 'load case: ', currentLoadCase - if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & - write(6,'(2x,a)') 'drop guessing along trajectory' - if (loadCases(currentLoadCase)%deformation%myType == 'l') then + if (.not. newLoadCase%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory' + if (newLoadCase%deformation%myType == 'l') then do j = 1_pInt, 3_pInt - if (any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & - any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .false.)) & - errorID = 832_pInt ! each row should be either fully or not at all defined + if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & + any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined enddo write(6,'(2x,a)') 'velocity gradient:' - else if (loadCases(currentLoadCase)%deformation%myType == 'f') then + else if (newLoadCase%deformation%myType == 'f') then write(6,'(2x,a)') 'deformation gradient at end of load case:' else write(6,'(2x,a)') 'deformation gradient rate:' endif do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt - if(loadCases(currentLoadCase)%deformation%maskLogical(i,j)) then - write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%deformation%values(i,j) + if(newLoadCase%deformation%maskLogical(i,j)) then + write(6,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j) else write(6,'(2x,12a)',advance='no') ' * ' endif enddo; write(6,'(/)',advance='no') enddo - if (any(loadCases(currentLoadCase)%stress%maskLogical .eqv. & - loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only - if (any(loadCases(currentLoadCase)%stress%maskLogical .and. & - transpose(loadCases(currentLoadCase)%stress%maskLogical) .and. & + if (any(newLoadCase%stress%maskLogical .eqv. & + newLoadCase%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only + if (any(newLoadCase%stress%maskLogical .and. & + transpose(newLoadCase%stress%maskLogical) .and. & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & errorID = 838_pInt ! no rotation is allowed by stress BC write(6,'(2x,a)') 'stress / GPa:' do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt - if(loadCases(currentLoadCase)%stress%maskLogical(i,j)) then - write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%stress%values(i,j)*1e-9_pReal + if(newLoadCase%stress%maskLogical(i,j)) then + write(6,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal else write(6,'(2x,12a)',advance='no') ' * ' endif enddo; write(6,'(/)',advance='no') enddo - if (any(abs(math_mul33x33(loadCases(currentLoadCase)%rotation, & - math_transpose33(loadCases(currentLoadCase)%rotation))-math_I3) > & + if (any(abs(math_mul33x33(newLoadCase%rotation, & + transpose(newLoadCase%rotation))-math_I3) > & reshape(spread(tol_math_check,1,9),[ 3,3]))& - .or. abs(math_det33(loadCases(currentLoadCase)%rotation)) > & + .or. abs(math_det33(newLoadCase%rotation)) > & 1.0_pReal + tol_math_check) errorID = 846_pInt ! given rotation matrix contains strain - if (any(dNeq(loadCases(currentLoadCase)%rotation, math_I3))) & + if (any(dNeq(newLoadCase%rotation, math_I3))) & write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& - math_transpose33(loadCases(currentLoadCase)%rotation) - if (loadCases(currentLoadCase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment - write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time - if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count - write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs - if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency - write(6,'(2x,a,i5)') 'output frequency: ', & - loadCases(currentLoadCase)%outputfrequency - write(6,'(2x,a,i5,/)') 'restart frequency: ', & - loadCases(currentLoadCase)%restartfrequency + transpose(newLoadCase%rotation) + if (newLoadCase%time < 0.0_pReal) errorID = 834_pInt ! negative time increment + write(6,'(2x,a,f12.6)') 'time: ', newLoadCase%time + if (newLoadCase%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count + write(6,'(2x,a,i5)') 'increments: ', newLoadCase%incs + if (newLoadCase%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency + write(6,'(2x,a,i5)') 'output frequency: ', newLoadCase%outputfrequency + write(6,'(2x,a,i5,/)') 'restart frequency: ', newLoadCase%restartfrequency if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message - enddo checkLoadcases - endif + endif reportAndCheck + loadCases = [loadCases,newLoadCase] ! load case is ok, append it + enddo + close(fileUnit) !-------------------------------------------------------------------------------------------------- -! doing initialization depending on selected solver +! doing initialization depending on active solvers call Utilities_init() do field = 1, nActiveFields select case (loadCases(1)%ID(field)) @@ -444,13 +432,13 @@ program DAMASK_spectral fileOffset = fileOffset + sum(outputSize) ! forward to current file position endif writeUndeformed !-------------------------------------------------------------------------------------------------- -! looping over loadcases +! looping over load cases loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) - time0 = time ! currentLoadCase start time + time0 = time ! load case start time guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc !-------------------------------------------------------------------------------------------------- -! loop over incs defined in input file for current currentLoadCase +! loop over incs defined in input file for current load case incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs totalIncsCounter = totalIncsCounter + 1_pInt @@ -460,13 +448,13 @@ program DAMASK_spectral if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) else - if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale - if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale + if (currentLoadCase == 1_pInt) then ! 1st load case of logarithmic scale + if (inc == 1_pInt) then ! 1st inc of 1st load case of logarithmic scale timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd - else ! not-1st inc of 1st currentLoadCase of logarithmic scale + else ! not-1st inc of 1st load case of logarithmic scale timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal)) endif - else ! not-1st currentLoadCase of logarithmic scale + else ! not-1st load case of logarithmic scale timeinc = time0 * & ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/& real(loadCases(currentLoadCase)%incs ,pReal))& @@ -633,8 +621,7 @@ program DAMASK_spectral convergedCounter, ' out of ', & notConvergedCounter + convergedCounter, ' (', & real(convergedCounter, pReal)/& - real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & - ' %) increments converged!' + real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!' flush(6) call MPI_file_close(resUnit,ierr) close(statUnit) @@ -655,10 +642,13 @@ end program DAMASK_spectral !-------------------------------------------------------------------------------------------------- subroutine quit(stop_id) #include - use MPI +#ifdef _OPENMP + use MPI, only: & + MPI_finalize +#endif use prec, only: & pInt - + implicit none integer(pInt), intent(in) :: stop_id integer, dimension(8) :: dateAndTime ! type default integer