From eeda357710526349fb94d715e09ad8e5e06c0797 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 30 Jan 2012 13:52:41 +0000 Subject: [PATCH] N^2 initialization loop (former spectralPictureMode) rewritten in material.f90 additional output in DAMASK_spectral_interface.f90 132 character cut off in constitutive_nonlocal.f90 rounding error in math.f90 complex number initialization (1.0_pReal)*2.0_pReal*pi new $DAMASK_NUM_THREADS warning in numerics.f90 / IO.f90 polishing in DAMASK_spectral.f90 --- code/DAMASK_spectral.f90 | 34 +++++++++------- code/DAMASK_spectral_interface.f90 | 28 +++++++++++--- code/IO.f90 | 2 + code/constitutive_nonlocal.f90 | 5 ++- code/material.f90 | 8 ++-- code/math.f90 | 2 +- code/numerics.f90 | 62 +++++++++++++++--------------- 7 files changed, 85 insertions(+), 56 deletions(-) diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index f8771df77..6927f8927 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -181,7 +181,8 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! variables for debugging fft using a scalar field - type(C_PTR) :: scalarField, plan_scalarField_forth, plan_scalarField_back + type(C_PTR) :: scalarField_realPointer, scalarField_complexPointer,& + plan_scalarField_forth, plan_scalarField_back real(pReal), dimension(:,:,:), pointer :: scalarField_real complex(pReal), dimension(:,:,:), pointer :: scalarField_complex integer(pInt) :: row, column @@ -199,7 +200,6 @@ program DAMASK_spectral print '(a,a)', ' Working Directory: ',trim(getSolverWorkingDirectoryName()) print '(a,a)', ' Solver Job Name: ',trim(getSolverJobName()) print '(a)', '' - !-------------------------------------------------------------------------------------------------- ! reading the load case file and allocate data structure containing load cases path = getLoadcaseName() @@ -387,7 +387,7 @@ program DAMASK_spectral print '(a,3(i12 ))','resolution a b c:', res print '(a,3(f12.5))','dimension x y z:', geomdim print '(a,i5)','homogenization: ',homog - if(any(cutting_freq/=0_pInt)) print '(a,3(i4),a)', 'cutting away ', cutting_freq, 'frequencies' + if(cut_off_value/=0.0_pReal) print '(a,3(i4),a)', 'cutting away ', cutting_freq, 'frequencies' print '(a)', '#############################################################' print '(a,a)', 'loadcase file: ',trim(getLoadcaseName()) @@ -557,12 +557,13 @@ program DAMASK_spectral endif if (debugFFTW) then - scalarField = fftw_alloc_complex(int(res1_red*res(2)*res(3),C_SIZE_T)) - call c_f_pointer(scalarField, scalarField_real, [ res(1)+2_pInt,res(2),res(3)]) - call c_f_pointer(scalarField, scalarField_complex, [ res1_red, res(2),res(3)]) - plan_scalarField_forth = fftw_plan_dft_r2c_3d(res(3),res(2),res(1),& !reversed order + scalarField_realPointer = fftw_alloc_complex(int(res(1) *res(2)*res(3),C_SIZE_T)) ! do not do an inplace transform + scalarField_complexPointer = fftw_alloc_complex(int(res1_red*res(2)*res(3),C_SIZE_T)) + call c_f_pointer(scalarField_realPointer, scalarField_real, [res(1), res(2),res(3)]) + call c_f_pointer(scalarField_complexPointer, scalarField_complex, [res1_red,res(2),res(3)]) + plan_scalarField_forth = fftw_plan_dft_r2c_3d(res(3),res(2),res(1),& !reversed order scalarField_real,scalarField_complex,fftw_planner_flag) - plan_scalarField_back = fftw_plan_dft_c2r_3d(res(3),res(2),res(1),& !reversed order + plan_scalarField_back = fftw_plan_dft_c2r_3d(res(3),res(2),res(1),& !reversed order scalarField_complex,scalarField_real,fftw_planner_flag) endif @@ -682,6 +683,7 @@ program DAMASK_spectral write(538), 'eoh' ! end of header write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results endif + if (debugGeneral) print '(a)' , 'Header of result file written out' if(totalIncsCounter > restartReadInc) then ! Do calculations (otherwise just forwarding) if(bc(loadcase)%restartFrequency>0_pInt) & @@ -739,7 +741,7 @@ program DAMASK_spectral c_reduced(k,j) = c_prev99(n,m) endif; enddo; endif; enddo call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) ! invert reduced stiffness - if(errmatinv) call IO_error(error_ID=800) + if(errmatinv) call IO_error(error_ID=799) s_prev99 = 0.0_pReal ! build full compliance k = 0_pInt do n = 1_pInt,9_pInt @@ -880,13 +882,19 @@ program DAMASK_spectral do k = 1_pInt, res(3); do j = 1_pInt, res(2) do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. err_div_RMS = err_div_RMS & - + 2.0_pReal*sum(abs(math_mul33x3_complex(tensorField_complex(i,j,k,1:3,1:3),& ! sum of squared absolute values of complex divergence. do not take square root and square again - xi(1:3,i,j,k))*differentationFactor)**2.0_pReal) ! --> sum squared L_2 norm of vector + + 2.0_pReal*(sum (real(math_mul33x3_complex(tensorField_complex(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again + xi(1:3,i,j,k))*differentationFactor)**2.0_pReal)&! --> sum squared L_2 norm of vector + +sum(aimag(math_mul33x3_complex(tensorField_complex(i,j,k,1:3,1:3),& + xi(1:3,i,j,k))*differentationFactor)**2.0_pReal)) enddo err_div_RMS = err_div_RMS & ! Those two layers do not have a conjugate complex counterpart - + sum(abs(math_mul33x3_complex(tensorField_complex(1 ,j,k,1:3,1:3),& + + sum(real(math_mul33x3_complex(tensorField_complex(1 ,j,k,1:3,1:3),& xi(1:3,1 ,j,k))*differentationFactor)**2.0_pReal)& - + sum(abs(math_mul33x3_complex(tensorField_complex(res1_red,j,k,1:3,1:3),& + + sum(aimag(math_mul33x3_complex(tensorField_complex(1 ,j,k,1:3,1:3),& + xi(1:3,1 ,j,k))*differentationFactor)**2.0_pReal)& + + sum(real(math_mul33x3_complex(tensorField_complex(res1_red,j,k,1:3,1:3),& + xi(1:3,res1_red,j,k))*differentationFactor)**2.0_pReal)& + + sum(aimag(math_mul33x3_complex(tensorField_complex(res1_red,j,k,1:3,1:3),& xi(1:3,res1_red,j,k))*differentationFactor)**2.0_pReal) enddo; enddo err_div_RMS = sqrt (err_div_RMS*wgt)* correctionFactor ! weighting by and taking square root (RMS). abs(...) because result is a complex number and multiplying with correction factor diff --git a/code/DAMASK_spectral_interface.f90 b/code/DAMASK_spectral_interface.f90 index 1602a6fb3..f289c6366 100644 --- a/code/DAMASK_spectral_interface.f90 +++ b/code/DAMASK_spectral_interface.f90 @@ -38,11 +38,11 @@ subroutine DAMASK_interface_init() implicit none - character(len=1024) commandLine + character(len=1024) commandLine, hostName, userName integer(pInt):: i, start = 0_pInt, length=0_pInt - + integer, dimension(8) :: date_and_time_values ! type default integer call get_command(commandLine) - + call DATE_AND_TIME(VALUES=date_and_time_values) do i=1,len(commandLine) ! remove capitals if(640) then ! if '--restart' is found, use that (contains '-l') + start = index(commandLine,'--restart',.true.) + 7_pInt + endif + length = index(commandLine(start:len(commandLine)),' ',.false.) + + call GET_ENVIRONMENT_VARIABLE('HOST',hostName) + call GET_ENVIRONMENT_VARIABLE('USER',userName) + write(6,*) write(6,*) '<<<+- DAMASK_spectral_interface init -+>>>' write(6,*) '$Id$' write(6,*) + write(6,'(a,2(i2.2,a),i4.4)'), ' Date: ',date_and_time_values(3),'/',& + date_and_time_values(2),'/',& + date_and_time_values(1) + write(6,'(a,2(i2.2,a),i2.2)'), ' Time: ',date_and_time_values(5),':',& + date_and_time_values(6),':',& + date_and_time_values(7) + write(6,*) 'Host Name: ', trim(hostName) + write(6,*) 'User Name: ', trim(userName) + write(6,*) 'Command line call: ', trim(commandLine) write(6,*) 'Geometry Parameter: ', trim(geometryParameter) write(6,*) 'Loadcase Parameter: ', trim(loadcaseParameter) - write(6,*) - !$OMP END CRITICAL (write2out) + if (start/=3_pInt) write(6,*) 'Restart Parameter: ', trim(commandLine(start:start+length)) endsubroutine DAMASK_interface_init diff --git a/code/IO.f90 b/code/IO.f90 index 752a69c4a..3c563545d 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1418,6 +1418,8 @@ endfunction select case (warning_ID) case (34) msg = 'invalid restart increment given' + case (35) + msg = 'could not get $DAMASK_NUM_THREADS' case (47_pInt) msg = 'No valid parameter for FFTW given, using FFTW_PATIENT' case (101_pInt) diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 2de72e4b9..c76bd450b 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -1859,8 +1859,9 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then endif if (considerLeavingFlux) then - normal_me2neighbor_defConf = math_det33(Favg) * math_mul33x3(math_inv33(math_transpose33(Favg)), mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) - normal_me2neighbor = math_mul33x3(math_transpose33(my_Fe), normal_me2neighbor_defConf) / math_det33(my_Fe) ! interface normal in my lattice configuration + normal_me2neighbor_defConf = math_det33(Favg) * math_mul33x3(math_inv33(math_transpose33(Favg)),& + mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) + normal_me2neighbor = math_mul33x3(math_transpose33(my_Fe), normal_me2neighbor_defConf) / math_det33(my_Fe) ! interface normal in my lattice configuration area = mesh_ipArea(n,ip,el) * math_norm3(normal_me2neighbor) normal_me2neighbor = normal_me2neighbor / math_norm3(normal_me2neighbor) ! normalize the surface normal to unit length do s = 1,ns diff --git a/code/material.f90 b/code/material.f90 index 79ac5dd0e..a4b5dc03d 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -234,7 +234,7 @@ subroutine material_parseHomogenization(file,myPart) allocate(homogenization_Noutput(Nsections)); homogenization_Noutput = 0_pInt allocate(homogenization_active(Nsections)); homogenization_active = .false. - forall (s = 1:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model? + forall (s = 1:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model? Homogenization view, maximum operations depend on maximum number of homog schemes homogenization_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) rewind(file) @@ -280,14 +280,14 @@ subroutine material_parseMicrostructure(file,myPart) use prec, only: pInt use IO - use mesh, only: mesh_element + use mesh, only: mesh_element, mesh_NcpElems implicit none character(len=*), intent(in) :: myPart integer(pInt), intent(in) :: file integer(pInt), parameter :: maxNchunks = 7 integer(pInt), dimension(1+2*maxNchunks) :: positions - integer(pInt) Nsections, section, constituent, i + integer(pInt) Nsections, section, constituent, e, i character(len=64) tag character(len=1024) line @@ -301,7 +301,7 @@ subroutine material_parseMicrostructure(file,myPart) allocate(microstructure_active(Nsections)) allocate(microstructure_elemhomo(Nsections)) - forall (i = 1:Nsections) microstructure_active(i) = any(mesh_element(4,:) == i) ! current microstructure used in model? + forall (e = 1:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements microstructure_Nconstituents = IO_countTagInPart(file,myPart,'(constituent)',Nsections) microstructure_maxNconstituents = maxval(microstructure_Nconstituents) diff --git a/code/math.f90 b/code/math.f90 index 27038e120..36bf948d6 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -3267,7 +3267,7 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) ! other variables integer(pInt) :: i, j, k, res1_red integer(pInt), dimension(3) :: k_s - complex(pReal), parameter :: integration_factor = cmplx(0.0_pReal,pi*2.0_pReal) + complex(pReal), parameter :: integration_factor = cmplx(0.0_pReal,1.0_pReal)*pi*2.0_pReal real(pReal), dimension(3) :: step, offset_coords if (debug_verbosity > 0_pInt) then diff --git a/code/numerics.f90 b/code/numerics.f90 index 5e5fe1c1a..21864815f 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -113,8 +113,9 @@ subroutine numerics_init() !*** output variables ***! !*** local variables ***! - integer(pInt), parameter :: fileunit = 300 - integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), parameter :: fileunit = 300_pInt + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt) :: gotDAMASK_NUM_THREADS = 1_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions character(len=64) tag character(len=1024) line @@ -135,48 +136,48 @@ subroutine numerics_init() iJacoStiffness = 1_pInt iJacoLpresiduum = 1_pInt pert_Fg = 1.0e-7_pReal - pert_method = 1 + pert_method = 1_pInt nHomog = 20_pInt subStepMinHomog = 1.0e-3_pReal - subStepSizeHomog = 0.25 - stepIncreaseHomog = 1.5 + subStepSizeHomog = 0.25_pReal + stepIncreaseHomog = 1.5_pReal nMPstate = 10_pInt nCryst = 20_pInt subStepMinCryst = 1.0e-3_pReal - subStepsizeCryst = 0.25 - stepIncreaseCryst = 1.5 + subStepsizeCryst = 0.25_pReal + stepIncreaseCryst = 1.5_pReal nState = 10_pInt nStress = 40_pInt rTol_crystalliteState = 1.0e-6_pReal rTol_crystalliteTemperature = 1.0e-6_pReal rTol_crystalliteStress = 1.0e-6_pReal - aTol_crystalliteStress = 1.0e-8_pReal ! residuum is in Lp (hence strain on the order of 1e-8 here) - numerics_integrator(1) = 1 ! fix-point iteration - numerics_integrator(2) = 1 ! fix-point iteration + aTol_crystalliteStress = 1.0e-8_pReal ! residuum is in Lp (hence strain on the order of 1e-8 here) + numerics_integrator(1) = 1_pInt ! fix-point iteration + numerics_integrator(2) = 1_pInt ! fix-point iteration !* RGC parameters: added <<>> with moderate setting - absTol_RGC = 1.0e+4 - relTol_RGC = 1.0e-3 - absMax_RGC = 1.0e+10 - relMax_RGC = 1.0e+2 - pPert_RGC = 1.0e-7 - xSmoo_RGC = 1.0e-5 - viscPower_RGC = 1.0e+0 ! Newton viscosity (linear model) - viscModus_RGC = 0.0e+0 ! No viscosity is applied - refRelaxRate_RGC = 1.0e-3 - maxdRelax_RGC = 1.0e+0 - maxVolDiscr_RGC = 1.0e-5 ! tolerance for volume discrepancy allowed - volDiscrMod_RGC = 1.0e+12 - volDiscrPow_RGC = 5.0 + absTol_RGC = 1.0e+4_pReal + relTol_RGC = 1.0e-3_pReal + absMax_RGC = 1.0e+10_pReal + relMax_RGC = 1.0e+2_pReal + pPert_RGC = 1.0e-7_pReal + xSmoo_RGC = 1.0e-5_pReal + viscPower_RGC = 1.0e+0_pReal ! Newton viscosity (linear model) + viscModus_RGC = 0.0e+0_pReal ! No viscosity is applied + refRelaxRate_RGC = 1.0e-3_pReal + maxdRelax_RGC = 1.0e+0_pReal + maxVolDiscr_RGC = 1.0e-5_pReal ! tolerance for volume discrepancy allowed + volDiscrMod_RGC = 1.0e+12_pReal + volDiscrPow_RGC = 5.0_pReal !* spectral parameters: - err_div_tol = 1.0e-4 ! 1.0e-4 proposed by Suquet - err_stress_tolrel = 0.01 ! relative tolerance for fullfillment of stress BC (1% of maximum stress) + err_div_tol = 1.0e-4_pReal ! 1.0e-4 proposed by Suquet + err_stress_tolrel = 0.01_pReal ! relative tolerance for fullfillment of stress BC (1% of maximum stress) itmax = 20_pInt ! Maximum iteration number memory_efficient = .true. ! Precalculate Gamma-operator (81 double per point) fftw_timelimit = -1.0_pReal ! no timelimit of plan creation for FFTW fftw_planner_string ='FFTW_PATIENT' - rotation_tol = 1.0e-12 + rotation_tol = 1.0e-12_pReal divergence_correction = .true. ! correct divergence by empirical factor simplified_algorithm = .true. ! use algorithm without fluctuation field update_gamma = .false. ! do not update gamma operator with current stiffness @@ -188,10 +189,11 @@ subroutine numerics_init() !* determin number of threads from environment variable DAMASK_NUM_THREADS DAMASK_NumThreadsInt = 0_pInt -!$ call GetEnv('DAMASK_NUM_THREADS',DAMASK_NumThreadsString) ! get environment variable DAMASK_NUM_THREADS... -!$ read(DAMASK_NumThreadsString,'(i4)') DAMASK_NumThreadsInt ! ...convert it to integer... -!$ if (DAMASK_NumThreadsInt < 1) DAMASK_NumThreadsInt = 1 ! ...ensure that its at least one... -!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! ...and use it as number of threads for parallel execution +!$ call GET_ENVIRONMENT_VARIABLE('DAMASK_NUM_THREADS',DAMASK_NumThreadsString,4_pInt,gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS... +!$ if(gotDAMASK_NUM_THREADS /= 0_pInt) call IO_warning(47,ext_msg=DAMASK_NumThreadsString) +!$ read(DAMASK_NumThreadsString,'(i4)') DAMASK_NumThreadsInt ! ...convert it to integer... +!$ if (DAMASK_NumThreadsInt < 1) DAMASK_NumThreadsInt = 1 ! ...ensure that its at least one... +!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! ...and use it as number of threads for parallel execution ! try to open the config file if(IO_open_file(fileunit,numerics_configFile)) then