223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

one more loop not needed

This commit is contained in:
Martin Diehl 2018-08-31 09:39:33 +02:00
parent 017563e061
commit 3cb279b083
1 changed files with 74 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

156
src/DAMASK_spectral.f90
View File

@ -115,7 +115,7 @@ program DAMASK_spectral
stagIterate stagIterate
integer(pInt) :: & integer(pInt) :: &
i, j, k, l, field, & i, j, k, l, field, &
errorID, & errorID = 0_pInt, &
cutBackLevel = 0_pInt, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ cutBackLevel = 0_pInt, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0_pInt !< fraction of current time interval stepFraction = 0_pInt !< fraction of current time interval
integer(pInt) :: & integer(pInt) :: &
@ -166,6 +166,7 @@ program DAMASK_spectral
if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1 if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1
if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1
allocate(solres(nActiveFields)) allocate(solres(nActiveFields))
allocate(newLoadCase%ID(nActiveFields))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assign mechanics solver depending on selected type ! assign mechanics solver depending on selected type
@ -188,7 +189,7 @@ program DAMASK_spectral
end select 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 allocate (loadCases(0)) ! array of load cases
open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read') 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)) 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 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 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 call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase
loadCases = [loadCases,newLoadCase]
currentLoadCase = currentLoadCase + 1_pInt newLoadCase%stress%myType='stress'
loadCases(currentLoadCase)%stress%myType='stress'
allocate(loadCases(currentLoadCase)%ID(nActiveFields))
field = 1 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 thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
field = field + 1 field = field + 1
loadCases(currentLoadCase)%ID(field) = FIELD_THERMAL_ID newLoadCase%ID(field) = FIELD_THERMAL_ID
endif thermalActive endif thermalActive
damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then
field = field + 1 field = field + 1
loadCases(currentLoadCase)%ID(field) = FIELD_DAMAGE_ID newLoadCase%ID(field) = FIELD_DAMAGE_ID
endif damageActive endif damageActive
do i = 1_pInt, chunkPos(1) do i = 1_pInt, chunkPos(1)
@ -230,46 +229,43 @@ program DAMASK_spectral
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'fdot'.or. & ! in case of Fdot, set type to fdot 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 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 else if (IO_lc(IO_stringValue(line,chunkPos,i)) == 'f') then
loadCases(currentLoadCase)%deformation%myType = 'f' newLoadCase%deformation%myType = 'f'
else else
loadCases(currentLoadCase)%deformation%myType = 'l' newLoadCase%deformation%myType = 'l'
endif endif
do j = 1_pInt, 9_pInt do j = 1_pInt, 9_pInt
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a * 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 if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo enddo
loadCases(currentLoadCase)%deformation%maskLogical = & ! logical mask in 3x3 notation newLoadCase%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) ! logical mask in 3x3 notation
transpose(reshape(temp_maskVector,[ 3,3])) newLoadCase%deformation%maskFloat = merge(ones,zeros,newLoadCase%deformation%maskLogical)! float (1.0/0.0) mask in 3x3 notation
loadCases(currentLoadCase)%deformation%maskFloat = & ! float (1.0/0.0) mask in 3x3 notation newLoadCase%deformation%values = math_plain9to33(temp_valueVector) ! values in 3x3 notation
merge(ones,zeros,loadCases(currentLoadCase)%deformation%maskLogical)
loadCases(currentLoadCase)%deformation%values = math_plain9to33(temp_valueVector) ! values in 3x3 notation
case('p','pk1','piolakirchhoff','stress', 's') case('p','pk1','piolakirchhoff','stress', 's')
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
do j = 1_pInt, 9_pInt do j = 1_pInt, 9_pInt
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk 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 if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo enddo
loadCases(currentLoadCase)%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) newLoadCase%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3]))
loadCases(currentLoadCase)%stress%maskFloat = merge(ones,zeros,& newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical)
loadCases(currentLoadCase)%stress%maskLogical) newLoadCase%stress%values = math_plain9to33(temp_valueVector)
loadCases(currentLoadCase)%stress%values = math_plain9to33(temp_valueVector)
case('t','time','delta') ! increment time 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 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) case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt) newLoadCase%incs = IO_intValue(line,chunkPos,i+1_pInt)
loadCases(currentLoadCase)%logscale = 1_pInt newLoadCase%logscale = 1_pInt
case('freq','frequency','outputfreq') ! frequency of result writings 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 case('r','restart','restartwrite') ! frequency of writing restart information
loadCases(currentLoadCase)%restartfrequency = & newLoadCase%restartfrequency = &
max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt)) max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt))
case('guessreset','dropguessing') case('guessreset','dropguessing')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory newLoadCase%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('euler') ! rotation of currentLoadCase given in euler angles case('euler') ! rotation of load case given in euler angles
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
l = 1_pInt ! assuming values given in degrees l = 1_pInt ! assuming values given in degrees
k = 1_pInt ! assuming keyword indicating degree/radians present 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) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j)
enddo enddo
if (l == 1_pInt) temp_valueVector(1:3) = temp_valueVector(1:3) * inRad ! convert to rad 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 newLoadCase%rotation = math_EulerToR(temp_valueVector(1:3)) ! convert rad Eulers to rotation matrix
case('rotation','rot') ! assign values for the rotation of currentLoadCase matrix case('rotation','rot') ! assign values for the rotation matrix
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
do j = 1_pInt, 9_pInt do j = 1_pInt, 9_pInt
temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j)
enddo enddo
loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) newLoadCase%rotation = math_plain9to33(temp_valueVector)
end select end select
enddo 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) reportAndCheck: if (worldrank == 0) then
!--------------------------------------------------------------------------------------------------
! 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)
write (loadcase_string, '(i6)' ) currentLoadCase write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & if (.not. newLoadCase%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory'
write(6,'(2x,a)') 'drop guessing along trajectory' if (newLoadCase%deformation%myType == 'l') then
if (loadCases(currentLoadCase)%deformation%myType == 'l') then
do j = 1_pInt, 3_pInt do j = 1_pInt, 3_pInt
if (any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. &
any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .false.)) & any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined
errorID = 832_pInt ! each row should be either fully or not at all defined
enddo enddo
write(6,'(2x,a)') 'velocity gradient:' 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:' write(6,'(2x,a)') 'deformation gradient at end of load case:'
else else
write(6,'(2x,a)') 'deformation gradient rate:' write(6,'(2x,a)') 'deformation gradient rate:'
endif endif
do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt
if(loadCases(currentLoadCase)%deformation%maskLogical(i,j)) then if(newLoadCase%deformation%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%deformation%values(i,j) write(6,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j)
else else
write(6,'(2x,12a)',advance='no') ' * ' write(6,'(2x,12a)',advance='no') ' * '
endif endif
enddo; write(6,'(/)',advance='no') enddo; write(6,'(/)',advance='no')
enddo enddo
if (any(loadCases(currentLoadCase)%stress%maskLogical .eqv. & if (any(newLoadCase%stress%maskLogical .eqv. &
loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only newLoadCase%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only
if (any(loadCases(currentLoadCase)%stress%maskLogical .and. & if (any(newLoadCase%stress%maskLogical .and. &
transpose(loadCases(currentLoadCase)%stress%maskLogical) .and. & transpose(newLoadCase%stress%maskLogical) .and. &
reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
errorID = 838_pInt ! no rotation is allowed by stress BC errorID = 838_pInt ! no rotation is allowed by stress BC
write(6,'(2x,a)') 'stress / GPa:' write(6,'(2x,a)') 'stress / GPa:'
do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt
if(loadCases(currentLoadCase)%stress%maskLogical(i,j)) then if(newLoadCase%stress%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%stress%values(i,j)*1e-9_pReal write(6,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal
else else
write(6,'(2x,12a)',advance='no') ' * ' write(6,'(2x,12a)',advance='no') ' * '
endif endif
enddo; write(6,'(/)',advance='no') enddo; write(6,'(/)',advance='no')
enddo enddo
if (any(abs(math_mul33x33(loadCases(currentLoadCase)%rotation, & if (any(abs(math_mul33x33(newLoadCase%rotation, &
math_transpose33(loadCases(currentLoadCase)%rotation))-math_I3) > & transpose(newLoadCase%rotation))-math_I3) > &
reshape(spread(tol_math_check,1,9),[ 3,3]))& 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 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:',& write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
math_transpose33(loadCases(currentLoadCase)%rotation) transpose(newLoadCase%rotation)
if (loadCases(currentLoadCase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment if (newLoadCase%time < 0.0_pReal) errorID = 834_pInt ! negative time increment
write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time write(6,'(2x,a,f12.6)') 'time: ', newLoadCase%time
if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count if (newLoadCase%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count
write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs write(6,'(2x,a,i5)') 'increments: ', newLoadCase%incs
if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency if (newLoadCase%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', & write(6,'(2x,a,i5)') 'output frequency: ', newLoadCase%outputfrequency
loadCases(currentLoadCase)%outputfrequency write(6,'(2x,a,i5,/)') 'restart frequency: ', newLoadCase%restartfrequency
write(6,'(2x,a,i5,/)') 'restart frequency: ', &
loadCases(currentLoadCase)%restartfrequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases endif reportAndCheck
endif 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() call Utilities_init()
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(1)%ID(field)) select case (loadCases(1)%ID(field))
@ -444,13 +432,13 @@ program DAMASK_spectral
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif writeUndeformed endif writeUndeformed
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! looping over loadcases ! looping over load cases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) 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 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 incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt totalIncsCounter = totalIncsCounter + 1_pInt
@ -460,13 +448,13 @@ program DAMASK_spectral
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else else
if (currentLoadCase == 1_pInt) then ! 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 currentLoadCase 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 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)) timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal))
endif endif
else ! not-1st currentLoadCase of logarithmic scale else ! not-1st load case of logarithmic scale
timeinc = time0 * & timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/& ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))& real(loadCases(currentLoadCase)%incs ,pReal))&
@ -633,8 +621,7 @@ program DAMASK_spectral
convergedCounter, ' out of ', & convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', & notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/& real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!'
' %) increments converged!'
flush(6) flush(6)
call MPI_file_close(resUnit,ierr) call MPI_file_close(resUnit,ierr)
close(statUnit) close(statUnit)
@ -655,7 +642,10 @@ end program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine quit(stop_id) subroutine quit(stop_id)
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
use MPI #ifdef _OPENMP
use MPI, only: &
MPI_finalize
#endif
use prec, only: & use prec, only: &
pInt pInt