From 5ebeb96e85daa92bd1bae2f4ab0a15ad341115fc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 6 Dec 2011 16:58:17 +0000 Subject: [PATCH] made DAMASK to work with gfortran: -removed to long lines -restructured f2py modules and merged make_DAMASK2Python into setup processing -setup_code.py now sets library path in makefile and asks for compile switches for spectral code -substituted \ in format strings with $ restructured DAMASK_spectral: -more logical output and structure of code -better input for spectral debug parameters --- {processing/post => code}/DAMASK.pyf | 0 code/DAMASK_spectral.f90 | 470 ++++++++++++----------- code/FEsolving.f90 | 32 +- code/IO.f90 | 4 +- code/constitutive_dislotwin.f90 | 6 +- code/crystallite.f90 | 3 +- code/debug.f90 | 17 +- code/makefile | 56 ++- code/math.f90 | 31 +- code/setup/setup_code.py | 26 +- processing/post/3Dvisualize.py | 4 +- processing/post/make_DAMASK2Python | 26 -- processing/pre/voronoi_randomSeeding.f90 | 4 +- processing/pre/voronoi_tessellation.f90 | 10 +- processing/setup/setup_processing.py | 137 +++++-- 15 files changed, 458 insertions(+), 368 deletions(-) rename {processing/post => code}/DAMASK.pyf (100%) delete mode 100644 processing/post/make_DAMASK2Python diff --git a/processing/post/DAMASK.pyf b/code/DAMASK.pyf similarity index 100% rename from processing/post/DAMASK.pyf rename to code/DAMASK.pyf diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 5113c82ae..8393d0618 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -52,7 +52,7 @@ program DAMASK_spectral use math use mesh, only: mesh_ipCenterOfGravity use CPFEM, only: CPFEM_general, CPFEM_initAll - use FEsolving, only: restartWrite, restartReadSpectral, restartReadStep + use FEsolving, only: restartWrite, restartReadStep use numerics, only: err_div_tol, err_stress_tolrel , rotation_tol,& itmax, memory_efficient, DAMASK_NumThreadsInt, divergence_correction, & fftw_planner_flag, fftw_timelimit @@ -103,10 +103,10 @@ program DAMASK_spectral logical :: spectralPictureMode = .false. ! indicating 1 to 1 mapping of FP to microstructure ! stress, stiffness and compliance average etc. - real(pReal), dimension(3,3) :: pstress, pstress_av, defgrad_av, & + real(pReal), dimension(3,3) :: pstress, pstress_av, defgrad_av_lab, & defgradAim = math_I3, defgradAimOld= math_I3, defgradAimCorr= math_I3,& mask_stress, mask_defgrad, fDot, & - pstress_av_load, defgradAim_lab ! quantities rotated to other coordinate system + pstress_av_lab, defgradAim_lab ! quantities rotated to other coordinate system real(pReal), dimension(3,3,3,3) :: dPdF, c0_reference, c_current = 0.0_pReal, s_prev, c_prev ! stiffness and compliance real(pReal), dimension(6) :: cstress ! cauchy stress real(pReal), dimension(6,6) :: dsde ! small strain stiffness @@ -149,7 +149,7 @@ program DAMASK_spectral bc_init%deformation = zeroes; bc_init%stress = zeroes; bc_init%rotation = zeroes bc_init%timeIncrement = 0.0_pReal; bc_init%temperature = 300.0_pReal bc_init%steps = 0_pInt; bc_init%logscale = 0_pInt - bc_init%outputfrequency = 1_pInt; bc_init%restartfrequency = 1_pInt + bc_init%outputfrequency = 1_pInt; bc_init%restartfrequency = 0_pInt bc_init%maskDeformation = .false.; bc_init%maskStress = .false. bc_init%maskStressVector = .false.; bc_init%velGradApplied = .false. bc_init%followFormerTrajectory = .true. @@ -242,7 +242,7 @@ program DAMASK_spectral case('f','freq','frequency','outputfreq') ! frequency of result writings bc(loadcase)%outputfrequency = IO_intValue(line,posLoadcase,j+1_pInt) case('r','restart','restartwrite') ! frequency of writing restart information - bc(loadcase)%restartfrequency = IO_intValue(line,posLoadcase,j+1_pInt) + bc(loadcase)%restartfrequency = max(0_pInt,IO_intValue(line,posLoadcase,j+1_pInt)) case('guessreset','dropguessing') bc(loadcase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory case('euler') ! rotation of loadcase given in euler angles @@ -367,27 +367,28 @@ program DAMASK_spectral if (.not. bc(loadcase)%followFormerTrajectory) print '(a)', 'drop guessing along trajectory' if (bc(loadcase)%velGradApplied) then do j = 1_pInt, 3_pInt - if (any(bc(loadcase)%maskDeformation(j,1:3) == .true.) .and. & - any(bc(loadcase)%maskDeformation(j,1:3) == .false.)) errorID = 32_pInt ! each row should be either fully or not at all defined + 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 enddo print '(a)','velocity gradient:' else print '(a)','deformation gradient rate:' endif - print '(3(3(f12.6,x)/)\)', merge(math_transpose3x3(bc(loadcase)%deformation),& + print '(3(3(f12.6,x)/)$)', merge(math_transpose3x3(bc(loadcase)%deformation),& reshape(spread(DAMASK_NaN,1,9),(/3,3/)),transpose(bc(loadcase)%maskDeformation)) - print '(a,/,3(3(f12.6,x)/)\)','stress / GPa:',1e-9*merge(math_transpose3x3(bc(loadcase)%stress),& + print '(a,/,3(3(f12.6,x)/)$)','stress / GPa:',1e-9*merge(math_transpose3x3(bc(loadcase)%stress),& reshape(spread(DAMASK_NaN,1,9),(/3,3/)),& transpose(bc(loadcase)%maskStress)) if (any(bc(loadcase)%rotation /= math_I3)) & - print '(a,3(3(f12.6,x)/)\)','rotation of loadframe:',math_transpose3x3(bc(loadcase)%rotation) + print '(a,3(3(f12.6,x)/)$)','rotation of loadframe:',math_transpose3x3(bc(loadcase)%rotation) print '(a,f12.6)','temperature:',bc(loadcase)%temperature print '(a,f12.6)','time: ',bc(loadcase)%timeIncrement print '(a,i5)','steps: ',bc(loadcase)%steps - print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency - print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency + print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency + print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency - if (any(bc(loadcase)%maskStress == bc(loadcase)%maskDeformation)) errorID = 31 ! exclusive or masking only + + if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 31 ! 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 @@ -398,8 +399,6 @@ program DAMASK_spectral if (bc(loadcase)%timeIncrement < 0.0_pReal) errorID = 34_pInt ! negative time increment if (bc(loadcase)%steps < 1_pInt) errorID = 35_pInt ! non-positive increment count if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 36_pInt ! non-positive result frequency - if (bc(loadcase)%restartfrequency < 1_pInt) errorID = 39_pInt ! non-positive restart frequency - if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) enddo @@ -455,9 +454,9 @@ program DAMASK_spectral !************************************************************* ! Initialization Start !************************************************************* - if(totalStepsCounter >= restartReadStep) then ! Do calculations (otherwise just forwarding) - - if (regrid==.true. ) then ! 'DeInitialize' the values changing in case of regridding + + if(totalStepsCounter == restartReadStep) then ! Initialize values + if (regrid .eqv. .true. ) then ! 'DeInitialize' the values changing in case of regridding regrid = .false. call dfftw_destroy_plan(fftw_plan(1)); call dfftw_destroy_plan(fftw_plan(2)) if(debugDivergence) call dfftw_destroy_plan(fftw_plan(3)) @@ -467,147 +466,153 @@ program DAMASK_spectral deallocate (temperature) deallocate (xi) deallocate (workfft) - !ToDo: here we have to create the new geometry and assign the values from the previous step + !ToDo: here we have to create the new geometry and assign the values from the previous step endif - if(totalStepsCounter == restartReadStep) then ! Initialize values - guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step - allocate (defgrad ( res(1),res(2),res(3),3,3)); defgrad = 0.0_pReal - allocate (defgradold ( res(1),res(2),res(3),3,3)); defgradold = 0.0_pReal - allocate (coordinates(3,res(1),res(2),res(3))); coordinates = 0.0_pReal - allocate (temperature( res(1),res(2),res(3))); temperature = bc(1)%temperature ! start out isothermally - allocate (xi (3,res(1)/2+1,res(2),res(3))); xi =0.0_pReal - allocate (workfft(res(1)+2,res(2),res(3),3,3)); workfft = 0.0_pReal - if (debugDivergence) allocate (divergence(res(1)+2,res(2),res(3),3)); divergence = 0.0_pReal + guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step + allocate (defgrad ( res(1),res(2),res(3),3,3)); defgrad = 0.0_pReal + allocate (defgradold ( res(1),res(2),res(3),3,3)); defgradold = 0.0_pReal + allocate (coordinates(3,res(1),res(2),res(3))); coordinates = 0.0_pReal + allocate (temperature( res(1),res(2),res(3))); temperature = bc(1)%temperature ! start out isothermally + allocate (xi (3,res(1)/2+1,res(2),res(3))); xi =0.0_pReal + allocate (workfft(res(1)+2,res(2),res(3),3,3)); workfft = 0.0_pReal + if (debugDivergence) allocate (divergence(res(1)+2,res(2),res(3),3)); divergence = 0.0_pReal - wgt = 1.0_pReal/real(Npoints, pReal) - call dfftw_plan_many_dft_r2c(fftw_plan(1),3,(/res(1),res(2),res(3)/),9,& - workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),& - workfft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),fftw_planner_flag) - call dfftw_plan_many_dft_c2r(fftw_plan(2),3,(/res(1),res(2),res(3)/),9,& - workfft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& - workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag) - if (debugDivergence) & - call dfftw_plan_many_dft_c2r(fftw_plan(3),3,(/res(1),res(2),res(3)/),3,& - divergence,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& - divergence,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag) - if (debugGeneral) then - !$OMP CRITICAL (write2out) - write (6,*) 'FFTW initialized' - !$OMP END CRITICAL (write2out) - endif - - if (restartReadStep==1_pInt) then ! not restarting, no deformation at the beginning - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - defgrad(i,j,k,1:3,1:3) = math_I3 - defgradold(i,j,k,1:3,1:3) = math_I3 - enddo; enddo; enddo - else ! using old values - if (IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(defgrad))) then - read (777,rec=1) defgrad - close (777) - endif - defgradold = defgrad - defgradAim = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - defgradAim = defgradAim + defgrad(i,j,k,1:3,1:3) ! calculating old average deformation - enddo; enddo; enddo - defgradAim = defgradAim * wgt - defgradAimOld = defgradAim - guessmode=0.0_pInt - endif - - ielem = 0_pInt - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - ielem = ielem + 1_pInt - coordinates(1:3,i,j,k) = mesh_ipCenterOfGravity(1:3,1,ielem) ! set to initial coordinates ToDo: SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!! But do we know them? I don't think so. Otherwise we don't need geometry reconstruction - call CPFEM_general(2_pInt,coordinates(1:3,i,j,k),math_I3,math_I3,temperature(i,j,k),0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) - c_current = c_current + dPdF - enddo; enddo; enddo - c0_reference = c_current * wgt ! linear reference material stiffness - c_prev = math_rotate_forward3x3x3x3(c0_reference,bc(loadcase)%rotation) ! rotate_forward: lab -> load system - - if (debugGeneral) then - !$OMP CRITICAL (write2out) - write (6,*) 'first call to CPFEM_general finished' - !$OMP END CRITICAL (write2out) - endif - - do k = 1_pInt, res(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - k_s(3) = k - 1_pInt - if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) - do j = 1_pInt, res(2) - k_s(2) = j - 1_pInt - if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) - do i = 1, res(1)/2_pInt + 1_pInt - k_s(1) = i - 1_pInt - xi(3,i,j,k) = 0.0_pReal ! 2D case - if(res(3) > 1_pInt) xi(3,i,j,k) = real(k_s(3), pReal)/geomdimension(3) ! 3D case - xi(2,i,j,k) = real(k_s(2), pReal)/geomdimension(2) - xi(1,i,j,k) = real(k_s(1), pReal)/geomdimension(1) - enddo; enddo; enddo - ! remove highest frequencies for calculation of divergence (CAREFUL, they will be used for pre calculatet gamma operator!) - do k = 1_pInt ,res(3); do j = 1_pInt ,res(2); do i = 1_pInt,res(1)/2_pInt + 1_pInt - if(k==res(3)/2_pInt+1_pInt) xi(3,i,j,k)= 0.0_pReal - if(j==res(2)/2_pInt+1_pInt) xi(2,i,j,k)= 0.0_pReal - if(i==res(1)/2_pInt+1_pInt) xi(1,i,j,k)= 0.0_pReal - enddo; enddo; enddo - - if(memory_efficient) then ! allocate just single fourth order tensor - allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal - else ! precalculation of gamma_hat field - allocate (gamma_hat(res(1)/2_pInt + 1_pInt ,res(2),res(3),3,3,3,3)); gamma_hat = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)/2_pInt + 1_pInt - if (any(xi(:,i,j,k) /= 0.0_pReal)) then - do l = 1_pInt ,3_pInt; do m = 1_pInt,3_pInt - xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k) - enddo; enddo - temp33_Real = math_inv3x3(math_mul3333xx33(c0_reference, xiDyad)) - else - xiDyad = 0.0_pReal - temp33_Real = 0.0_pReal - endif - do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt; do p=1_pInt,3_pInt - gamma_hat(i,j,k, l,m,n,p) = - 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *& - (xiDyad(m,p)+xiDyad(p,m)) - enddo; enddo; enddo; enddo - enddo; enddo; enddo - endif - - ! write header of output file + wgt = 1.0_pReal/real(Npoints, pReal) + call dfftw_plan_many_dft_r2c(fftw_plan(1),3,(/res(1),res(2),res(3)/),9,& + workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),& + workfft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),fftw_planner_flag) + call dfftw_plan_many_dft_c2r(fftw_plan(2),3,(/res(1),res(2),res(3)/),9,& + workfft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& + workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag) + if (debugDivergence) & + call dfftw_plan_many_dft_c2r(fftw_plan(3),3,(/res(1),res(2),res(3)/),3,& + divergence,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& + divergence,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag) + if (debugGeneral) then !$OMP CRITICAL (write2out) - open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())& - //'.spectralOut',form='UNFORMATTED',status='REPLACE') - write(538), 'load', trim(getLoadcaseName()) - write(538), 'workingdir', trim(getSolverWorkingDirectoryName()) - write(538), 'geometry', trim(getSolverJobName())//InputFileExtension - write(538), 'resolution', res - write(538), 'dimension', geomdimension - write(538), 'materialpoint_sizeResults', materialpoint_sizeResults - write(538), 'loadcases', N_Loadcases - write(538), 'logscale', bc(1:N_Loadcases)%logscale ! one entry per loadcase (0: linear, 1: log) - write(538), 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase - write(538), 'times', bc(1:N_Loadcases)%timeIncrement ! one entry per loadcase - bc(1)%steps= bc(1)%steps + 1_pInt - write(538), 'increments', bc(1:N_Loadcases)%steps ! one entry per loadcase ToDo: rename keyword to steps - bc(1)%steps= bc(1)%steps - 1_pInt - write(538), 'startingIncrement', restartReadStep -1_pInt ! start with writing out the previous step - write(538), 'eoh' ! end of header - write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed) results - !$OMP END CRITICAL (write2out) + write (6,*) 'FFTW initialized' + !$OMP END CRITICAL (write2out) endif + + if (restartReadStep==1_pInt) then ! not restarting, no deformation at the beginning + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + defgrad(i,j,k,1:3,1:3) = math_I3 + defgradold(i,j,k,1:3,1:3) = math_I3 + enddo; enddo; enddo + else ! using old values + if (IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(defgrad))) then + read (777,rec=1) defgrad + close (777) + endif + defgradold = defgrad + defgradAim = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + defgradAim = defgradAim + defgrad(i,j,k,1:3,1:3) ! calculating old average deformation + enddo; enddo; enddo + defgradAim = defgradAim * wgt + defgradAimOld = defgradAim + guessmode=0.0_pInt + endif + ielem = 0_pInt + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ielem = ielem + 1_pInt + coordinates(1:3,i,j,k) = mesh_ipCenterOfGravity(1:3,1,ielem) ! set to initial coordinates ToDo: SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!! But do we know them? I don't think so. Otherwise we don't need geometry reconstruction + call CPFEM_general(2_pInt,coordinates(1:3,i,j,k),math_I3,math_I3,temperature(i,j,k),& + 0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) + c_current = c_current + dPdF + enddo; enddo; enddo + c0_reference = c_current * wgt ! linear reference material stiffness + + if (debugGeneral) then + !$OMP CRITICAL (write2out) + write (6,*) 'first call to CPFEM_general finished' + !$OMP END CRITICAL (write2out) + endif + + do k = 1_pInt, res(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) + k_s(3) = k - 1_pInt + if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) + do j = 1_pInt, res(2) + k_s(2) = j - 1_pInt + if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) + do i = 1, res(1)/2_pInt + 1_pInt + k_s(1) = i - 1_pInt + xi(3,i,j,k) = 0.0_pReal ! 2D case + if(res(3) > 1_pInt) xi(3,i,j,k) = real(k_s(3), pReal)/geomdimension(3) ! 3D case + xi(2,i,j,k) = real(k_s(2), pReal)/geomdimension(2) + xi(1,i,j,k) = real(k_s(1), pReal)/geomdimension(1) + enddo; enddo; enddo + !remove highest frequencies for calculation of divergence (CAREFULL, they will be used for pre calculatet gamma operator!) + do k = 1_pInt ,res(3); do j = 1_pInt ,res(2); do i = 1_pInt,res(1)/2_pInt + 1_pInt + if(k==res(3)/2_pInt+1_pInt) xi(3,i,j,k)= 0.0_pReal + if(j==res(2)/2_pInt+1_pInt) xi(2,i,j,k)= 0.0_pReal + if(i==res(1)/2_pInt+1_pInt) xi(1,i,j,k)= 0.0_pReal + enddo; enddo; enddo + + if(memory_efficient) then ! allocate just single fourth order tensor + allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal + else ! precalculation of gamma_hat field + allocate (gamma_hat(res(1)/2_pInt + 1_pInt ,res(2),res(3),3,3,3,3)); gamma_hat = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)/2_pInt + 1_pInt + if (any(xi(:,i,j,k) /= 0.0_pReal)) then + do l = 1_pInt ,3_pInt; do m = 1_pInt,3_pInt + xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k) + enddo; enddo + temp33_Real = math_inv3x3(math_mul3333xx33(c0_reference, xiDyad)) + else + xiDyad = 0.0_pReal + temp33_Real = 0.0_pReal + endif + do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt; do p=1_pInt,3_pInt + gamma_hat(i,j,k, l,m,n,p) = - 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *& + (xiDyad(m,p)+xiDyad(p,m)) + enddo; enddo; enddo; enddo + enddo; enddo; enddo + endif + + if (divergence_correction) then + if (res(3) == 1_pInt) then + correctionFactor = minval(geomdimension(1:2))*wgt**(-1.0_pReal/4.0_pReal) ! 2D case, ToDo: correct? + else + correctionFactor = minval(geomdimension(1:3))*wgt**(-1.0_pReal/4.0_pReal) ! multiplying by minimum dimension to get rid of dimension dependency and phenomenologigal factor wgt**(-1/4) to get rid of resolution dependency + endif + else + correctionFactor = 1.0_pReal + endif + + ! write header of output file + !$OMP CRITICAL (write2out) + open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())& + //'.spectralOut',form='UNFORMATTED',status='REPLACE') + write(538), 'load', trim(getLoadcaseName()) + write(538), 'workingdir', trim(getSolverWorkingDirectoryName()) + write(538), 'geometry', trim(getSolverJobName())//InputFileExtension + write(538), 'resolution', res + write(538), 'dimension', geomdimension + write(538), 'materialpoint_sizeResults', materialpoint_sizeResults + write(538), 'loadcases', N_Loadcases + write(538), 'logscale', bc(1:N_Loadcases)%logscale ! one entry per loadcase (0: linear, 1: log) + write(538), 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase + write(538), 'times', bc(1:N_Loadcases)%timeIncrement ! one entry per loadcase + bc(1)%steps= bc(1)%steps + 1_pInt + write(538), 'increments', bc(1:N_Loadcases)%steps ! one entry per loadcase ToDo: rename keyword to steps + bc(1)%steps= bc(1)%steps - 1_pInt + write(538), 'startingIncrement', restartReadStep -1_pInt ! start with writing out the previous step + write(538), 'eoh' ! end of header + write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed) results + !$OMP END CRITICAL (write2out) + endif !************************************************************* ! Initialization End !************************************************************* - - if (mod(step - 1_pInt,bc(loadcase)%restartFrequency)==0_pInt) then ! at frequency of writing restart information - restartWrite = .true. ! setting restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?) - else - restartWrite = .false. - endif - - if (bc(loadcase)%velGradApplied) & ! calculate fDot from given L and current F - fDot = math_mul33x33(bc(loadcase)%deformation, defgradAim) + + if(totalStepsCounter >= restartReadStep) then ! Do calculations (otherwise just forwarding) + if(bc(loadcase)%restartFrequency>0_pInt) & + restartWrite = ( mod(step - 1_pInt,bc(loadcase)%restartFrequency)==0_pInt) ! at frequency of writing restart information + ! setting restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?) + if (bc(loadcase)%velGradApplied) & ! calculate fDot from given L and current F + fDot = math_mul33x33(bc(loadcase)%deformation, defgradAim) !winding forward of deformation aim in loadcase system temp33_Real = defgradAim @@ -645,7 +650,8 @@ program DAMASK_spectral CPFEM_mode = 1_pInt ! winding forward iter = 0_pInt err_div = 2.0_pReal * err_div_tol ! go into loop - + + c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness from former step if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied c_prev99 = math_Plain3333to99(c_prev) k = 0_pInt ! build reduced stiffness @@ -674,17 +680,17 @@ program DAMASK_spectral s_prev = (math_Plain99to3333(s_prev99)) endif - !$OMP CRITICAL (write2out) + print '(a)', '#############################################################' print '(A,I5.5,A,es12.6)', 'Increment ', totalStepsCounter, ' Time ',time if (restartWrite ) then print '(A)', 'writing converged results of previous step for restart' - if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then ! and writing deformation gradient field to file - write (777,rec=1) defgrad - close (777) - endif - endif - !$OMP END CRITICAL (write2out) + if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then ! and writing deformation gradient field to file + write (777,rec=1) defgrad + close (777) + endif + endif + !************************************************************* ! convergence loop do while(iter < itmax .and. & @@ -692,34 +698,27 @@ program DAMASK_spectral err_stress > err_stress_tol)) iter = iter + 1_pInt !************************************************************* + print '(a)', '' print '(a)', '=============================================================' print '(5(a,i5.5))', 'Loadcase ',loadcase,' Step ',step,'/',bc(loadcase)%steps,'@Iteration ',iter,'/',itmax do n = 1_pInt,3_pInt; do m = 1_pInt,3_pInt - defgrad_av(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n)) * wgt + defgrad_av_lab(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n)) * wgt enddo; enddo - !$OMP CRITICAL (write2out) - print '(a,/,3(3(f12.7,x)/)\)', 'deformation gradient:',math_transpose3x3(defgrad_av) - print '(l)', restartWrite + print '(a,/,3(3(f12.7,x)/)$)', 'deformation gradient:',& + math_transpose3x3(math_rotate_forward3x3(defgrad_av_lab,bc(loadcase)%rotation)) + print '(a)', '' print '(a)', '... update stress field P(F) ................................' - !$OMP END CRITICAL (write2out) - defgradDetMax = -huge(1.0_pReal) - defgradDetMin = +huge(1.0_pReal) + ielem = 0_pInt do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - defgradDet = math_det3x3(defgrad(i,j,k,1:3,1:3)) - defgradDetMax = max(defgradDetMax,defgradDet) - defgradDetMin = min(defgradDetMin,defgradDet) ielem = ielem + 1_pInt - call CPFEM_general(3_pInt,& ! collect cycle + call CPFEM_general(3_pInt,& ! collect cycle coordinates(1:3,i,j,k), defgradold(i,j,k,1:3,1:3), defgrad(i,j,k,1:3,1:3),& temperature(i,j,k),timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) enddo; enddo; enddo - print '(a,x,es10.4)' , 'max determinant of deformation:', defgradDetMax - print '(a,x,es10.4)' , 'min determinant of deformation:', defgradDetMin - workfft = 0.0_pReal ! needed because of the padding for FFTW c_current = 0.0_pReal ielem = 0_pInt @@ -731,43 +730,24 @@ program DAMASK_spectral temperature(i,j,k),timeinc,ielem,1_pInt,& cstress,dsde, pstress,dPdF) CPFEM_mode = 2_pInt - workfft(i,j,k,1:3,1:3) = pstress ! build up average P-K stress + + workfft(i,j,k,1:3,1:3) = pstress ! build up average P-K stress c_current = c_current + dPdF enddo; enddo; enddo - restartWrite = .false. ! ToDo: don't know if we need it. Depends on how CPFEM_general is writing results + do n = 1_pInt,3_pInt; do m = 1_pInt,3_pInt - pstress_av(m,n) = sum(workfft(1:res(1),1:res(2),1:res(3),m,n)) * wgt + pstress_av_lab(m,n) = sum(workfft(1:res(1),1:res(2),1:res(3),m,n)) * wgt enddo; enddo - !$OMP CRITICAL (write2out) - print '(a,/,3(3(f12.7,x)/)\)', 'Piola-Kirchhoff stress / MPa: ',math_transpose3x3(pstress_av)/1.e6 - - err_stress_tol = 0.0_pReal - pstress_av_load = math_rotate_forward3x3(pstress_av,bc(loadcase)%rotation) - - if(size_reduced > 0_pInt) then ! calculate stress BC if applied - err_stress = maxval(abs(mask_stress * (pstress_av_load - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable) - err_stress_tol = maxval(abs(mask_defgrad * pstress_av_load)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent - print '(a)', '' - print '(a,es10.4,a,f6.2)', 'error stress = ',err_stress, ', ', err_stress/err_stress_tol - print '(a)', '... correcting deformation gradient to fulfill BCs ..........' - defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av_load - bc(loadcase)%stress))) ! residual on given stress components - defgradAim = defgradAim + defgradAimCorr - print '(a,/,3(3(f12.7,x)/)\)', 'new deformation aim: ',& - math_transpose3x3(math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation)) - print '(a,x,es10.4)' , 'with determinant: ', math_det3x3(defgradAim) - endif - print '(a)', '' print '(a)', '... calculating equilibrium with spectral method ............' - !$OMP END CRITICAL (write2out) call dfftw_execute_dft_r2c(fftw_plan(1),workfft,workfft) ! FFT of pstress p_hat_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(workfft(1,1,1,1:3,1:3),& ! L_2 norm of average stress in fourier space, math_transpose3x3(workfft(1,1,1,1:3,1:3)))))) ! ignore imaginary part as it is always zero for real only input)) err_div = 0.0_pReal - err_div_max = 0.0_pReal + err_div_max = 0.0_pReal ! only if if(debugDivergence) == .true. of importace do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)/2_pInt+1_pInt err_div = err_div + sqrt(sum((& ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain) math_mul33x3_complex(workfft(i*2_pInt-1_pInt,j,k,1:3,1:3) + & @@ -781,17 +761,9 @@ program DAMASK_spectral xi(1:3,i,j,k))& )**2.0_pReal)))) enddo; enddo; enddo - if (divergence_correction) then - if (res(3) == 1_pInt) then - correctionFactor = minval(geomdimension(1:2))*wgt**(-1.0_pReal/4.0_pReal) ! 2D case, ToDo: correct? PE: Do we need this in the loop or can be pre-calculated? - else - correctionFactor = minval(geomdimension(1:3))*wgt**(-1.0_pReal/4.0_pReal) ! multiplying by minimum dimension to get rid of dimension dependency and phenomenologigal factor wgt**(-1/4) to get rid of resolution dependency - endif - else - correctionFactor = 1.0_pReal - endif + err_div = err_div*wgt/p_hat_avg*correctionFactor ! weighting by points and average stress and multiplying with correction factor - err_div_max = err_div_max/p_hat_avg*correctionFactor ! weighting by average stress and multiplying with correction factor + err_div_max = err_div_max/p_hat_avg*correctionFactor ! weighting by average stress and multiplying with correction factor, only if if(debugDivergence) == .true. of importace if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res(1)/2_pInt+1_pInt @@ -825,6 +797,8 @@ program DAMASK_spectral workfft(i*2_pInt ,j,k,1:3,1:3) = aimag(temp33_Complex) enddo; enddo; enddo endif + + ! calculate additional divergence criteria and report ------------- if(debugDivergence) then divergence = 0.0_pReal do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)/2_pInt+1_pInt @@ -845,27 +819,11 @@ program DAMASK_spectral err_real_div_avg = err_real_div_avg + sqrt(sum((divergence(i,j,k,1:3))**2.0_pReal)) ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain) err_real_div_max = max(err_real_div_max,abs(sqrt(sum((divergence(i,j,k,1:3))**2.0_pReal)))) ! maximum of L two norm of div(stress) in fourier space (Suquet large strain) enddo; enddo; enddo - p_real_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(pstress_av,& ! L_2 norm of average stress in fourier space, - math_transpose3x3(pstress_av))))) ! ignore imaginary part as it is always zero for real only input)) + p_real_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(pstress_av_lab,& ! L_2 norm of average stress in real space, + math_transpose3x3(pstress_av_lab))))) ! ignore imaginary part as it is always zero for real only input)) err_real_div_avg = err_real_div_avg*wgt/p_real_avg err_real_div_max = err_real_div_max/p_real_avg - endif - ! average strain - workfft(1,1,1,1:3,1:3) = defgrad_av - math_I3 ! zero frequency (real part) - workfft(2,1,1,1:3,1:3) = 0.0_pReal ! zero frequency (imaginary part) - - call dfftw_execute_dft_c2r(fftw_plan(2),workfft,workfft) - defgrad = defgrad + workfft(1:res(1),:,:,:,:)*wgt - do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt - defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt - enddo; enddo - defgradAim_lab = math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation) - do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt - defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim_lab(m,n) - defgrad_av(m,n)) ! anticipated target minus current state - enddo; enddo - !$OMP CRITICAL (write2out) - if(debugDivergence) then print '(a,es10.4,a,f6.2)', 'error divergence FT avg = ',err_div, ', ', err_div/err_div_tol print '(a,es10.4)', 'error divergence FT max = ',err_div_max print '(a,es10.4)', 'error divergence Real avg = ',err_real_div_avg @@ -873,12 +831,63 @@ program DAMASK_spectral else print '(a,es10.4,a,f6.2)', 'error divergence = ',err_div, ', ', err_div/err_div_tol endif - !$OMP END CRITICAL (write2out) + ! -------------------------- + + workfft(1,1,1,1:3,1:3) = defgrad_av_lab - math_I3 ! assign zero frequency (real part) with average displacement gradient + workfft(2,1,1,1:3,1:3) = 0.0_pReal ! zero frequency (imaginary part) + + call dfftw_execute_dft_c2r(fftw_plan(2),workfft,workfft) ! backtransformation of fluct deformation gradient + + defgrad = defgrad + workfft(1:res(1),1:res(2),1:res(3),1:3,1:3)*wgt ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization + + do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt + defgrad_av_lab(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n))*wgt ! ToDo: check whether this needs recalculation or is equivalent to former defgrad_av + enddo; enddo + + + ! stress boundary condition check ------------- + pstress_av = math_rotate_forward3x3(pstress_av_lab,bc(loadcase)%rotation) + print '(a,/,3(3(f12.7,x)/)$)', 'Piola-Kirchhoff stress / MPa: ',math_transpose3x3(pstress_av)/1.e6 + + if(size_reduced > 0_pInt) then ! calculate stress BC if applied + err_stress = maxval(abs(mask_stress * (pstress_av - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable) + err_stress_tol = maxval(abs(mask_defgrad * pstress_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent + print '(a)', '' + print '(a,es10.4,a,f6.2)', 'error stress = ',err_stress, ', ', err_stress/err_stress_tol + print '(a)', '... correcting deformation gradient to fulfill BCs ..........' + defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av - bc(loadcase)%stress))) ! residual on given stress components + defgradAim = defgradAim + defgradAimCorr + print '(a,/,3(3(f12.7,x)/)$)', 'new deformation aim: ', math_transpose3x3(defgradAim) + print '(a,x,es10.4)' , 'with determinant: ', math_det3x3(defgradAim) + else + err_stress_tol = 0.0_pReal + endif + ! ------------------------------ + + ! homogeneous correction towards avg deformation gradient ------------- + defgradAim_lab = math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation) ! boundary conditions from load frame into lab (Fourier) frame + do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt + defgrad(1:res(1),1:res(2),1:res(3),m,n) = & + defgrad(1:res(1),1:res(2),1:res(3),m,n) + (defgradAim_lab(m,n) - defgrad_av_lab(m,n)) ! anticipated target minus current state + enddo; enddo + ! ------------------------------ + + ! bounds of det(F) ------------- + defgradDetMax = -huge(1.0_pReal) + defgradDetMin = +huge(1.0_pReal) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + defgradDet = math_det3x3(defgrad(i,j,k,1:3,1:3)) + defgradDetMax = max(defgradDetMax,defgradDet) + defgradDetMin = min(defgradDetMin,defgradDet) + enddo; enddo; enddo + + print '(a,x,es10.4)' , 'max determinant of deformation:', defgradDetMax + print '(a,x,es10.4)' , 'min determinant of deformation:', defgradDetMin + + ! ------------------------------ enddo ! end looping when convergency is achieved - c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness for next step - !ToDo: Incfluence for next loadcase !$OMP CRITICAL (write2out) print '(a)', '' print '(a)', '=============================================================' @@ -894,7 +903,7 @@ program DAMASK_spectral write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file endif !$OMP END CRITICAL (write2out) - endif + endif ! end calculation/forwarding enddo ! end looping over steps in current loadcase deallocate(c_reduced) deallocate(s_reduced) @@ -902,7 +911,8 @@ program DAMASK_spectral !$OMP CRITICAL (write2out) print '(a)', '' print '(a)', '#############################################################' - print '(a,i5.5,a,i5.5,a)', 'of ', totalStepsCounter - restartReadStep + 1_pInt, ' calculated steps, ', notConvergedCounter, ' steps did not converge!' + print '(a,i5.5,a,i5.5,a)', 'of ', totalStepsCounter - restartReadStep + 1_pInt, ' calculated steps, ',& + notConvergedCounter, ' steps did not converge!' !$OMP END CRITICAL (write2out) close(538) call dfftw_destroy_plan(fftw_plan(1)); call dfftw_destroy_plan(fftw_plan(2)) diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index 1dede9637..ea0fb484b 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -62,28 +62,22 @@ FEmodelGeometry = getModelName() if (IO_open_inputFile(fileunit,FEmodelGeometry)) then if (trim(FEsolver) == 'Spectral') then - call get_command(commandLine) ! may contain uppercase + 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 + commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase enddo - start = index(commandLine,'-r ',.true.) + 3_pInt ! search for '-r' and jump forward to given name - if (index(commandLine,'--restart ',.true.)>0) then ! if '--restart' is found, use that (contains '-r') - start = index(commandLine,'--restart ',.true.) + 10_pInt - endif - length = index(commandLine(start:len(commandLine)),' ',.false.) - if(start /= 3_pInt) then ! found at least -r - read(commandLine(start:start+length),'(I12)') restartReadStep + 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)') restartReadStep ! read argument + restartRead = max(0_pInt,restartReadStep) > 0_pInt + if(restartReadStep < 0_pInt) call IO_warning(warning_ID=34_pInt) endif - if (restartReadStep > 0_pInt) then - restartRead = .true. - restartWrite = .true. - endif - if (restartReadStep == 0_pInt) then - restartRead = .false. - restartWrite = .false. - endif - if(restartReadStep < 0_pInt) call IO_error(error_ID=47) else rewind(fileunit) do @@ -142,7 +136,7 @@ endif enddo elseif (FEsolver == 'Spectral') then - restartReadSpectral = .true. + !do nothing else call IO_error(106) ! cannot open file for old job info endif diff --git a/code/IO.f90 b/code/IO.f90 index 2259ad2a5..07ac4eb63 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1166,8 +1166,6 @@ endfunction msg = 'incomplete information in spectral mesh header' case (46) msg = 'not a rotation defined for loadcase rotation' - case (47) - msg = 'invalid restart increment given' case (50) msg = 'writing constitutive output description' case (100) @@ -1416,6 +1414,8 @@ endfunction select case (warning_ID) case (33_pInt) msg = 'cannot guess along trajectory for first step of first loadcase' + case (34) + msg = 'invalid restart increment given' case (47_pInt) msg = 'No valid parameter for FFTW given, using FFTW_PATIENT' case (101_pInt) diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index ca52dd28c..0ff685ba0 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -1057,7 +1057,8 @@ if(constitutive_dislotwin_sbVelocity(myInstance) /= 0.0_pReal) then !* Stress ratios StressRatio_p = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(myInstance))**constitutive_dislotwin_p(myInstance) - StressRatio_pminus1 = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(myInstance))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) + StressRatio_pminus1 = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(myInstance))& + **(constitutive_dislotwin_p(myInstance)-1.0_pReal) !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) !* Initial shear rates @@ -1463,7 +1464,8 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) tau = dot_product(Tstar_v,constitutive_dislotwin_sbSv(1:6,j,g,ip,el)) !* Stress ratios StressRatio_p = (abs(tau)/constitutive_dislotwin_sbResistance(myInstance))**constitutive_dislotwin_p(myInstance) - StressRatio_pminus1 = (abs(tau)/constitutive_dislotwin_sbResistance(myInstance))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) + StressRatio_pminus1 = (abs(tau)/constitutive_dislotwin_sbResistance(myInstance))& + **(constitutive_dislotwin_p(myInstance)-1.0_pReal) !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) !* Initial shear rates diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 61a537af8..4c31c7353 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3065,7 +3065,8 @@ logical error ! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL --- ! --- we use crystallite_orientation from above, so need a seperate loop -!$OMP PARALLEL DO PRIVATE(myPhase,myInstance,myStructure,neighboring_e,neighboring_i,neighboringPhase,neighboringInstance,neighboringStructure) +!$OMP PARALLEL DO PRIVATE(myPhase,myInstance,myStructure,neighboring_e,neighboring_i,neighboringPhase,& +!$OMP neighboringInstance,neighboringStructure) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) myPhase = material_phase(1,i,e) ! get my phase diff --git a/code/debug.f90 b/code/debug.f90 index 550ffbce3..6e1cb4757 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -115,7 +115,7 @@ subroutine debug_init() select case(tag) case ('element','e','el') debug_e = IO_intValue(line,positions,2) - case ('itegrationpoint','i','ip') + case ('integrationpoint','i','ip') debug_i = IO_intValue(line,positions,2) case ('grain','g','gr') debug_g = IO_intValue(line,positions,2) @@ -123,12 +123,15 @@ subroutine debug_init() debug_selectiveDebugger = IO_intValue(line,positions,2) > 0_pInt case ('verbosity') debug_verbosity = IO_intValue(line,positions,2) - case ('(generaldebugspectral)') ! use bitwise logical and, continue with +8_pInt - if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 1_pInt - case ('(divergencedebugspectral)') - if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 2_pInt - case ('(restartdebugspectral)') - if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 4_pInt + case ('(spectral)') + select case(IO_lc(IO_stringValue(line,positions,2))) + case('general') + spectral_debug_verbosity = ior(spectral_debug_verbosity, 1) + case('divergence') + spectral_debug_verbosity = ior(spectral_debug_verbosity, 2) + case('restart') + spectral_debug_verbosity = ior(spectral_debug_verbosity, 4) + endselect endselect enddo 100 close(fileunit) diff --git a/code/makefile b/code/makefile index 1be6d486f..2d78f5def 100644 --- a/code/makefile +++ b/code/makefile @@ -25,10 +25,13 @@ # PORTABLE = TRUE (FALSE): decision, if executable is optimized for the machine on which it was built. # OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE,ULTRA): Optimization mode: O2, O0, O3 + further options for most files, 03 + further options for all files # OPENMP = TRUE (FALSE): OpenMP multiprocessor support -# ACMLPATH =/opt/acml4.4.0/ifort64/lib (...) Path to ACML Library, choose according to your system -# ACMLPATH =/opt/acml4.4.0/ifort64_mp/lib (...) Path to ACML Library with multicore support, choose according to your system -# FFTWOPTIONS =./../lib/libfftw3.a include/libfftw3_threads.a -lpthread (...) Path to FFTW library files with Linux threads (multicore) support -# FFTWOPTIONS =./../lib/libfftw3.a (...) Path to FFTW library files without Linux threads (multicore) support +# FFTWPATH =TAKE_FFTW_PATH, will be adjusted by setup_code.py +# ACMLROOT =TAKE_ACLM_ROOT, will be adjusted by setup_code.py +# ACMLPATH =$(ACMLROOT)/"compilerdir"/lib (...) Path to ACML Library, choose according to your system +# ACMLPATH =$(ACMLROOT/"compilerdir"_mp/lib (...) Path to ACML Library with multicore support, choose according to your system +# "compilerdir" is "intel64" for ifort +# FFTWOPTIONS =$(FFTWPATH)/libfftw3.a $(FFTWPATH)/libfftw3_threads.a -lpthread (...) Path to FFTW library files with Linux threads (multicore) support +# FFTWOPTIONS =$(FFTWPATH)/libfftw3.a (...) Path to FFTW library files without Linux threads (multicore) support # FFTWOPTIONS is different for single and double precision. Choose the options to use OpenMP instead of pthreads support or change the directory # PREFIX: specify an arbitrary prefix # SUFFIX: specify an arbitrary suffix @@ -38,17 +41,29 @@ # information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/ # check if an array index is too small (<1) or too large! DEBUG1 =-check bounds -g + #will cause a lot of warnings because we create a bunch of temporary arrays DEBUG2 =-check arg_temp_created + #check from time to time DEBUG3 =-fp-stack-check -g -traceback -gen-interfaces -warn interfaces + #should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Problably it helps also to unlimit other limits DEBUG4 =-heap-arrays + #checks for standard DEBUG5 =stand std03/std95 -#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3) + +#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3) ######################################################################################## +#default values below will be set by setup_code.py +#FFTWPATH =../lib +#ACMLROOT =/opt/acml4.4.0 + +FFTWPATH =/nethome/m.diehl/DAMASK/lib +ACMLROOT =/opt/acml4.4.0 + ifndef F90 F90 =ifort endif @@ -75,34 +90,35 @@ ifndef OPENMP OPENMP=ON endif +# setting defaults in case of multicore support ifeq ($(OPENMP),ON) OPENMP_FLAG_ifort =-openmp -openmp-report0 -parallel OPENMP_FLAG_gfortran =-fopenmp ifndef ACMLPATH -ACMLPATH =/opt/acml4.4.0/ifort64_mp/lib +ACMLPATH =$(ACMLROOT)/$(F90)64_mp/lib endif ifndef FFTWOPTIONS ifeq ($(PRECISION),single) -FFTWOPTIONS =./../lib/libfftw3f_threads.a ./../lib/libfftw3f.a -lpthread +FFTWOPTIONS =$(FFTWPATH)/libfftw3f_threads.a $(FFTWPATH)/libfftw3f.a -lpthread else -FFTWOPTIONS =./../lib/libfftw3_threads.a ./../lib/libfftw3.a -lpthread +FFTWOPTIONS =$(FFTWPATH)/libfftw3_threads.a $(FFTWPATH)/libfftw3.a -lpthread endif endif -BLAS_ifort =-L $(ACMLPATH) -lacml_mp -BLAS_gfortran = +BLAS=-L $(ACMLPATH) -lacml_mp + +#setting defaults in case of single core compilation else ifndef ACMLPATH -ACMLPATH =/opt/acml4.4.0/ifort64/lib +ACMLPATH=$(ACMLROOT)/$(F90)64/lib endif ifndef FFTWOPTIONS ifeq ($(PRECISION),single) -FFTWOPTIONS =./../lib/libfftw3f.a +FFTWOPTIONS =$(FFTWPATH)/libfftw3f.a else -FFTWOPTIONS =./../lib/libfftw3.a +FFTWOPTIONS =$(FFTWPATH)/libfftw3.a endif endif -BLAS_ifort =-L $(ACMLPATH) -lacml -BLAS_gfortran = +BLAS=-L $(ACMLPATH) -lacml endif @@ -111,10 +127,10 @@ OPTIMIZATION_OFF_gfortran =-O0 OPTIMIZATION_DEFENSIVE_ifort =-O2 OPTIMIZATION_DEFENSIVE_gfortran =-O2 OPTIMIZATION_AGGRESSIVE_ifort =-O3 $(PORTABLE_SWITCH) -ip -static -fp-model fast=2 -no-prec-div -OPTIMIZATION_AGGRESSIVE_gfortran =-O3 $(PORTABLE_SWITCH) -march=opteron -ffast-math -funroll-loops -ftree-vectorize -ftree-loop-linear +OPTIMIZATION_AGGRESSIVE_gfortran =-O3 $(PORTABLE_SWITCH) -march=opteron -ffast-math -funroll-loops -ftree-vectorize -ftree-loop-linear COMPILE_OPTIONS_ifort =-fpp -diag-disable 8291,8290 -COMPILE_OPTIONS_gfortran =-xf95-cpp-input -ffree-line-length-none -fno-range-check +COMPILE_OPTIONS_gfortran =-xf95-cpp-input -fno-range-check COMPILE =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${OPTIMIZATION}_${F90}} -c COMPILE_MAXOPTI =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${MAXOPTI}_${F90}} -c @@ -125,14 +141,14 @@ endif ifeq ($(PRECISION),single) DAMASK_spectral_single.exe: DAMASK_spectral_single.o CPFEM.a - $(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral_single.exe DAMASK_spectral_single.o CPFEM.a $(FFTWOPTIONS)\ - constitutive.a advanced.a basics.a ${BLAS_${F90}} + $(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) -o DAMASK_spectral_single.exe DAMASK_spectral_single.o CPFEM.a $(FFTWOPTIONS)\ + constitutive.a advanced.a basics.a $(BLAS) DAMASK_spectral_single.o: DAMASK_spectral_single.f90 CPFEM.o $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) DAMASK_spectral_single.f90 $(SUFFIX) else DAMASK_spectral.exe: DAMASK_spectral.o CPFEM.a $(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a $(FFTWOPTIONS)\ - constitutive.a advanced.a basics.a ${BLAS_${F90}} + constitutive.a advanced.a basics.a $(BLAS) DAMASK_spectral.o: DAMASK_spectral.f90 CPFEM.o $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) DAMASK_spectral.f90 $(SUFFIX) endif diff --git a/code/math.f90 b/code/math.f90 index e20b1eadc..2503199a4 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -3048,7 +3048,8 @@ subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes) me = (/i,j,k/) ! me on skin shift = sign(abs(res+diag-2_pInt*me)/(res+diag),res+diag-2_pInt*me) lookup = me-diag+shift*res - wrappedCentroids(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = centroids(lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt,1:3) - & + wrappedCentroids(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = & + centroids(lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt,1:3) - & matmul(defgrad_av, shift*geomdim) endif enddo; enddo; enddo @@ -3056,7 +3057,8 @@ subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes) do j = 0_pInt,res(2) do i = 0_pInt,res(1) do n = 1_pInt,8_pInt - nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = nodes(i+1_pInt,j+1_pInt,k+1_pInt,3) + wrappedCentroids(i+1_pInt+neighbor(1_pInt,n), & + nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = & + nodes(i+1_pInt,j+1_pInt,k+1_pInt,3) + wrappedCentroids(i+1_pInt+neighbor(1_pInt,n), & j+1_pInt+neighbor(2,n), & k+1_pInt+neighbor(3,n),1:3) enddo; enddo; enddo; enddo @@ -3197,6 +3199,7 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) ! other variables integer(pInt) :: i, j, k integer(pInt), dimension(3) :: k_s + complex(pReal), parameter :: integration_factor = cmplx(0.0_pReal,pi*2.0_pReal) real(pReal), dimension(3) :: step, offset_coords integer*8, dimension(2) :: plan_fft @@ -3228,11 +3231,11 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) do i = 1_pInt, res(1)/2_pInt+1_pInt k_s(1) = i-1_pInt if(i/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)& - + defgrad_fft(i,j,k,1:3,1)*geomdim(1)/(real(k_s(1),pReal)*cmplx(0.0_pReal,1.0_pReal)*pi*2.0_pReal) + + defgrad_fft(i,j,k,1:3,1)*geomdim(1)/(real(k_s(1),pReal)*integration_factor) if(j/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)& - + defgrad_fft(i,j,k,1:3,2)*geomdim(2)/(real(k_s(2),pReal)*cmplx(0.0_pReal,1.0_pReal)*pi*2.0_pReal) + + defgrad_fft(i,j,k,1:3,2)*geomdim(2)/(real(k_s(2),pReal)*integration_factor) if(k/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)& - + defgrad_fft(i,j,k,1:3,3)*geomdim(3)/(real(k_s(3),pReal)*cmplx(0.0_pReal,1.0_pReal)*pi*2.0_pReal) + + defgrad_fft(i,j,k,1:3,3)*geomdim(3)/(real(k_s(3),pReal)*integration_factor) enddo; enddo; enddo call dfftw_execute_dft_c2r(plan_fft(2), coords_fft, coords) @@ -3432,12 +3435,18 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence_field) order = order + 1_pInt do k = 0_pInt, res(3)-1_pInt; do j = 0_pInt, res(2)-1_pInt; do i = 0_pInt, res(1)-1_pInt do m = 1_pInt, order - coordinates(1,1:3) = mesh_location(mesh_index((/i+m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) - coordinates(2,1:3) = mesh_location(mesh_index((/i-m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) - coordinates(3,1:3) = mesh_location(mesh_index((/i,j+m,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) - coordinates(4,1:3) = mesh_location(mesh_index((/i,j-m,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) - coordinates(5,1:3) = mesh_location(mesh_index((/i,j,k+m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) - coordinates(6,1:3) = mesh_location(mesh_index((/i,j,k-m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) + coordinates(1,1:3) = mesh_location(mesh_index((/i+m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& + + (/1_pInt,1_pInt,1_pInt/) + coordinates(2,1:3) = mesh_location(mesh_index((/i-m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& + + (/1_pInt,1_pInt,1_pInt/) + coordinates(3,1:3) = mesh_location(mesh_index((/i,j+m,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& + + (/1_pInt,1_pInt,1_pInt/) + coordinates(4,1:3) = mesh_location(mesh_index((/i,j-m,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& + + (/1_pInt,1_pInt,1_pInt/) + coordinates(5,1:3) = mesh_location(mesh_index((/i,j,k+m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& + + (/1_pInt,1_pInt,1_pInt/) + coordinates(6,1:3) = mesh_location(mesh_index((/i,j,k-m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& + + (/1_pInt,1_pInt,1_pInt/) do l = 1_pInt, vec_tens divergence_field(i+1_pInt,j+1_pInt,k+1_pInt,l) = divergence_field(i+1_pInt,j+1_pInt,k+1_pInt,l) + FDcoefficient(m,order) * & ((field(coordinates(1,1),coordinates(1,2),coordinates(1,3),l,1)- & diff --git a/code/setup/setup_code.py b/code/setup/setup_code.py index 01c34aad7..a62073914 100755 --- a/code/setup/setup_code.py +++ b/code/setup/setup_code.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # $Id$ # Writes version specific files for different MARC releases -import os,sys +import os,sys,string,damask_tools architectures = { 'marc': { @@ -29,8 +29,28 @@ for arch in architectures: childFile.write(line.replace(me['versions'][0],version)) childFile.close() -if raw_input("Do you want to compile the spectral code now? (y/n)") is 'y': - os.system('make --directory %s'%wd) +# changing dirs in make file +damask_variables = damask_tools.DAMASK_TOOLS() +makefile = open(os.path.join(damask_variables.rootDir(),'code/makefile')) +content = makefile.readlines() +makefile.close() +makefile = open(os.path.join(damask_variables.rootDir(),'code/makefile'),'w') +for line in content: + if line.startswith('FFTWPATH'): + line='FFTWPATH =%s\n'%(damask_variables.pathInfo['fftw']) + print line + if line.startswith('ACMLROOT'): + line='ACMLROOT =%s\n'%(damask_variables.pathInfo['acml']) + print line + makefile.writelines(line) +makefile.close() + +# compiling spectral code +if raw_input("Do you want to compile the spectral code now? (y/n) ") is 'y' or 'Y': + compiler_switches = raw_input("Please give compiling switches (Enter to use default) ") + os.system('make --directory %s clean'%(wd)) + print compiler_switches + os.system('make --directory %s %s'%(wd,compiler_switches)) if '--clean' in [s.lower() for s in sys.argv]: os.system('make --directory %s clean'%wd) diff --git a/processing/post/3Dvisualize.py b/processing/post/3Dvisualize.py index 458e688e6..a4292e84c 100755 --- a/processing/post/3Dvisualize.py +++ b/processing/post/3Dvisualize.py @@ -310,7 +310,7 @@ parser.add_option('--box', dest='output_box', action='store_true', \ help='produce VTK box file') parser.add_option('--nobox', dest='output_box', action='store_false', \ help='omit VTK box file') -parser.add_option('--scaling', dest='scaling', type='int', \ +parser.add_option('--scaling', dest='scaling', type='float', \ help='scaling of fluctuation [%default]') parser.add_option('-u', '--unitlength', dest='unitlength', type='float', \ help='set unit length for 2D model [%default]') @@ -321,7 +321,7 @@ parser.set_defaults(tensor = []) parser.set_defaults(output_mesh = True) parser.set_defaults(output_points = False) parser.set_defaults(output_box = False) -parser.set_defaults(scaling = 1) +parser.set_defaults(scaling = 1.0) parser.set_defaults(undeformed = False) parser.set_defaults(unitlength = 0.0) parser.set_defaults(cell = True) diff --git a/processing/post/make_DAMASK2Python b/processing/post/make_DAMASK2Python deleted file mode 100644 index 26fac525d..000000000 --- a/processing/post/make_DAMASK2Python +++ /dev/null @@ -1,26 +0,0 @@ -#!/bin/bash - -# $Id: make_postprocessingMath 1106 2011-11-17 21:36:56Z MPIE\m.diehl $ - -# This script is used to compile math.f90 and make the functions defined in DAMASK_math.pyf -# avialable for python in the module DAMASK_math.so -# It uses the fortran wrapper f2py that is included in the numpy package to construct the -# module postprocessingMath.so out of the fortran code postprocessingMath.f90 -# for the generation of the pyf file: -#f2py -m DAMASK -h DAMASK.pyf --overwrite-signature ../../code/math.f90 \ - -if [[ $# -eq 0 ]]; then - wd='.' -else - wd=$1 -fi -# use env. variables here! -cd $wd -rm ../../lib/DAMASK.so -f2py -c \ -DAMASK.pyf \ -../../code/DAMASK2Python_helper.f90 \ -../../code/math.f90 \ -../../lib/libfftw3.a \ -/opt/acml4.4.0/ifort64/lib/libacml.a -mv DAMASK.so ../../lib/. diff --git a/processing/pre/voronoi_randomSeeding.f90 b/processing/pre/voronoi_randomSeeding.f90 index a5d467a72..4e754eb37 100644 --- a/processing/pre/voronoi_randomSeeding.f90 +++ b/processing/pre/voronoi_randomSeeding.f90 @@ -81,8 +81,8 @@ program voronoi open(21, file = trim(filename)//('.seeds')) write(21, '(i1,a1,a6)') 4,achar(9),'header' write(21, '(A, I8, A, I8, A, I8)') 'resolution a ', a, ' b ', b, ' c ', c - write(21,*) 'grains', N_Seeds - write(21,*) 'random seed ',rndInit(1) + write(21, '(A, I8)') 'grains', N_Seeds + write(21, '(A, I8)') 'random seed ',rndInit(1) write(21,'(6(a,a1))') 'x',achar(9),'y',achar(9),'z',achar(9),'phi1',achar(9),'Phi',achar(9),'phi2',achar(9) do i = 1, n_Seeds write(21, '(6(F10.6,a1))'),seeds(i,1), achar(9), seeds(i,2), achar(9), seeds(i,3), achar(9), & diff --git a/processing/pre/voronoi_tessellation.f90 b/processing/pre/voronoi_tessellation.f90 index 52c1d69b0..7aaf7c7d8 100644 --- a/processing/pre/voronoi_tessellation.f90 +++ b/processing/pre/voronoi_tessellation.f90 @@ -169,7 +169,7 @@ program voronoi use IO implicit none - logical gotN_Seeds, gotResolution + logical :: gotN_Seeds=.false., gotResolution = .false. logical, dimension(:), allocatable :: grainCheck character(len=1024) input_name, output_name, format1, format2, N_Digits, line, key integer(pInt) N_Seeds, theGrain, i, j, k, l, m @@ -273,14 +273,14 @@ program voronoi if(geomdim(i) .gt. 0.0) validDim(i) = .true. enddo - if(all(validDim == .false.)) then + if(all(validDim .eqv. .false.)) then geomdim(maxval(maxloc(resolution))) = 1.0 validDim(maxval(maxloc(resolution))) = .true. print*, 'no valid dimension specified, using automated setting' endif do i=1,3 - if (validDim(i) == .false.) then + if (validDim(i) .eqv. .false.) then print*, 'rescaling ivalid dimension' , i geomdim(i) = maxval(geomdim/real(resolution),validDim)*real(resolution(i)) endif @@ -317,11 +317,11 @@ program voronoi open(21, file = ((trim(output_name))//'.spectral')) write(20, '(A)'), '3 header' write(20, '(A, I8, A, I8, A, I8)'), 'resolution a ', resolution(1), ' b ', resolution(2), ' c ', resolution(3) - write(20, '(A, g15.10, A, g15.10, A, g15.10)'), 'dimension x ', geomdim(1), ' y ', geomdim(2), ' z ', geomdim(3) + write(20, '(A, g17.10, A, g17.10, A, g17.10)'), 'dimension x ', geomdim(1), ' y ', geomdim(2), ' z ', geomdim(3) write(20, '(A)'), 'homogenization 1' format1 = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' ! geom format - format2 = '(3(tr2, f6.2), 3(tr2,g10.5), I10, a)' ! spectral (Lebensohn) format + format2 = '(3(tr2, f6.2), 3(tr2,g10.5), I10, a)' ! spectral (Lebensohn) format ! perform voronoi tessellation and write result to files diff --git a/processing/setup/setup_processing.py b/processing/setup/setup_processing.py index 293fe2133..febdf262b 100755 --- a/processing/setup/setup_processing.py +++ b/processing/setup/setup_processing.py @@ -2,11 +2,62 @@ # Makes postprocessing routines acessible from everywhere. -import os,sys,glob +import os,sys,glob,string +from optparse import OptionParser, Option -import damask_tools; reload(damask_tools) -damask_tools.DAMASK_TOOLS().check_env() +import damask_tools +# ----------------------------- +class extendableOption(Option): +# ----------------------------- +# used for definition of new option parser action 'extend', which enables to take multiple option arguments +# taken from online tutorial http://docs.python.org/library/optparse.html + + ACTIONS = Option.ACTIONS + ("extend",) + STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",) + TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",) + ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",) + + def take_action(self, action, dest, opt, value, values, parser): + if action == "extend": + lvalue = value.split(",") + values.ensure_value(dest, []).extend(lvalue) + else: + Option.take_action(self, action, dest, opt, value, values, parser) + + + +######################################################## +# MAIN +######################################################## + +parser = OptionParser(option_class=extendableOption, usage='%prog options', description = """ +Sets up the pre and post processing tools of DAMASK + +""" + string.replace('$Id: addDivergence.py 1129 2011-12-01 12:01:13Z MPIE\m.diehl $','\n','\\n') +) + +parser.add_option('--F90', dest='compiler', type='string', \ + help='name of F90 compiler [%default]') + +parser.set_defaults(compiler = 'ifort') +(options,filenames) = parser.parse_args() + +#translating name of compiler for use with f2py and setting subdirname of acml +if options.compiler == 'gfortran': + f2py_compiler='gnu95 --f90flags="-fno-range-check"' + acml_subdir='ifort64/lib' +else: + f2py_compiler='intelem' + acml_subdir='ifort64/lib' + + +#getting pathinfo +damask_variables = damask_tools.DAMASK_TOOLS() +baseDir = os.path.join(damask_variables.rootDir(),'processing/') +codeDir = os.path.join(damask_variables.rootDir(),'code/') + +#define ToDo list bin_link = { \ 'pre' : [ 'marc_addUserOutput.py', @@ -16,8 +67,8 @@ bin_link = { \ 'spectral_geomCheck.py', 'spectral_minimalSurface.py', 'spectral_vicinityOffset.py', - 'voronoi_randomSeeding.py', - 'voronoi_tessellation.py', + 'voronoi_randomSeeding.exe', + 'voronoi_tessellation.exe', ], 'post' : [ '3Dvisualize.py', @@ -46,28 +97,58 @@ compile = { \ 'post' : [ ], } + execute = { \ - 'pre' : [ - ], - 'post' : [ - 'make_DAMASK2Python', - ], + 'postMath' : [ + 'rm %slib/DAMASK.so' %(os.path.join(damask_variables.rootDir(),'lib/')), + # The following command is used to compile math.f90 and make the functions defined in DAMASK_math.pyf + # avialable for python in the module DAMASK_math.so + # It uses the fortran wrapper f2py that is included in the numpy package to construct the + # module postprocessingMath.so out of the fortran code postprocessingMath.f90 + # for the generation of the pyf file: + #f2py -m DAMASK -h DAMASK.pyf --overwrite-signature ../../code/math.f90 \ + 'f2py %sDAMASK.pyf '%(codeDir) +\ + '-c --fcompiler=%s '%(f2py_compiler) +\ + '%sDAMASK2Python_helper.f90 ' %(codeDir)+\ + '%smath.f90 ' %(codeDir)+\ + '%s/libfftw3.a ' %(damask_variables.pathInfo['fftw'])+\ + '%s/%s/libacml.a' %(damask_variables.pathInfo['acml'],acml_subdir), + 'mv %sDAMASK.so %s/.' %(codeDir,os.path.join(damask_variables.rootDir(),'lib/')), + ] } -damask_root = os.getenv('DAMASK_ROOT') -if damask_root == '' or damask_root == None: damask_root = os.path.join(os.path.dirname(sys.argv[0]),'../../') -damask_bin = os.getenv('DAMASK_BIN') -if damask_bin == '' or damask_bin == None: damask_bin = os.path.join(damask_root,'bin/') -basedir = os.path.join(damask_root,'processing/') +for dir in compile: + for file in compile[dir]: + src = os.path.abspath(os.path.join(baseDir,dir,file)) + if os.path.isfile(src): + try: + os.system('rm %s.exe'%(os.path.splitext(src)[0])) + print 'removing %s.exe '%(os.path.splitext(src)[0]) + except: + pass + print 'compiling ',src,'using',options.compiler + os.system('%s -O2 -o%s.exe %s'%(options.compiler,os.path.splitext(src)[0],src)) + +os.chdir(codeDir) # needed for compilation with gfortran and f2py +for tasks in execute: + for cmd in execute[tasks]: + os.system(cmd) +os.chdir(os.path.join(damask_variables.rootDir(),'processing/setup/')) + +modules = glob.glob('*.mod') +for module in modules: + print 'removing', module + os.remove(module) + for dir in bin_link: for file in bin_link[dir]: - src = os.path.abspath(os.path.join(basedir,dir,file)) + src = os.path.abspath(os.path.join(baseDir,dir,file)) if (file == ''): - sym_link = os.path.abspath(os.path.join(damask_bin,dir)) + sym_link = os.path.abspath(os.path.join(damask_variables.binDir(),dir)) else: - sym_link = os.path.abspath(os.path.join(damask_bin,os.path.splitext(file)[0])) + sym_link = os.path.abspath(os.path.join(damask_variables.binDir(),os.path.splitext(file)[0])) print sym_link,'-->',src if os.path.lexists(sym_link): os.remove(sym_link) @@ -78,23 +159,3 @@ for dir in bin_link: #if os.path.lexists(old_link): # os.remove(old_link) - -for dir in compile: - for file in compile[dir]: - src = os.path.abspath(os.path.join(basedir,dir,file)) - if os.path.isfile(src): - print file - os.system('rm %s.exe'%(os.path.splitext(src)[0])) - os.system('ifort -O3 -parallel -o%s %s'%(os.path.splitext(src)[0],src)) - -modules = glob.glob('*.mod') -for module in modules: - print module - os.remove(module) - -for dir in execute: - for file in execute[dir]: - cmd = os.path.abspath(os.path.join(basedir,dir,file)) - if os.access(cmd,os.X_OK): - print cmd - os.system('%s %s'%(cmd,os.path.join(basedir,dir)))