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
This commit is contained in:
Martin Diehl
parent
7731926fec
commit
eeda357710
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32)
|
||||
enddo
|
||||
|
|
@ -117,15 +117,31 @@ subroutine DAMASK_interface_init()
|
|||
loadcaseParameter = '' ! should be empty
|
||||
loadcaseParameter(1:length)=commandLine(start:start+length)
|
||||
|
||||
!$OMP CRITICAL (write2out)
|
||||
start = index(commandLine,'-r',.true.) + 3_pInt ! search for '-r' and jump forward iby 3 to given name
|
||||
if (index(commandLine,'--restart',.true.)>0) 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
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 <<<updated 17.12.2009>>> 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
|
||||
|
|
|
|||
Loading…
Reference in New Issue